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1  INTRODUCTION 


Nascap-2k  is  a  spacecraft  charging  and  plasma  interactions  code  designed  to  be  used  by 
spacecraft  designers,  aerospace  and  materials  engineers,  and  space  plasma  environments  experts 
to  study  the  effects  of  both  the  natural  and  spacecraft-generated  plasma  environments  on 
spacecraft  systems. 

Designers  of  spacecraft  for  government,  commercial,  and  research  purposes  require  advanced 
modeling  capabilities  to  guide  the  design  of  satellites  that  can  survive  and  operate  properly  in  the 
natural  environment.  Computer  modeling  of  flight  experiments  (such  as  SCATHA  (Spacecraft 
Charging  at  High  Altitude),  [1]  the  SPEAR  (Space  Power  Experiments  Aboard  Rockets)  [2,  3] 
series,  and  CHAWS  (Charging  Hazards  and  Wake  Studies) [4])  has  demonstrated  excellent  ability 
to  predict  both  steady- state  and  dynamic  interactions  between  high-voltage  spacecraft  and  the 
ambient  plasma.  This  ability  was  also  extended  to  inherently  dynamic  problems  involving  three- 
dimensional  space  charge  sheath  formation,  current  flow  in  the  quasi-neutral  presheath, 
breakdown  phenomena,  plasma  kinetics,  ionization  processes,  and  the  effect  of  unsteady 
processes  on  spacecraft  charging. 

Nascap-2k  encapsulates  the  knowledge  gained  in  these  efforts  as  well  as  in  modeling  spacecraft 
with  unique  requirements  and  destined  for  disparate  environments.  This  gives  the  spacecraft 
designer  quality  modeling  capabilities  by  taking  advantage  of  the  present  understanding  of  the 
pertinent  phenomena,  employing  advanced  algorithms,  and  implementing  a  state-of-the-art  user 
interface,  including  three-dimensional  post-processing  graphics. 

The  core  capabilities  of  Nascap-2k  are  as  follows: 

■  Define  spacecraft  surfaces  and  geometry  and  the  structure  of  the  computational  space 
surrounding  the  spacecraft. 

■  Solve  for  time-dependent  potentials  on  spacecraft  surfaces. 

■  Solve  the  electrostatic  potential  around  the  object,  with  flexible  boundary  conditions  on  the 
object  and  with  space-charge  computed  either  fully  by  particles,  fully  analytically,  or  in  a 
hybrid  manner. 

■  Generate,  track,  and  otherwise  process  particles  of  various  species,  represented  as  macro¬ 
particles,  in  the  computational  space. 

■  View  surface  potentials,  space  potentials,  particle  trajectories,  and  time -dependent  potentials 
and  currents. 

Nascap-2k  calculates  surface  charging  in  tenuous  plasma  environments,  such  as  geosynchronous 
earth  orbit  (GEO)  and  interplanetary  (Solar  Wind)  environments,  using  the  Boundary  Element 
Method  [5]  (BEM).  The  Boundary  Element  Method  facilitates  calculation  of  surface  electric 
fields  (which  limit  the  emission  of  photoelectrons  and  secondary  electrons)  without  the  need  to 
grid  the  space  surrounding  the  spacecraft.  It  also  enables  Nascap-2k  to  anticipate  electric  field 
changes  due  to  surface  charging,  resulting  in  a  smoother  and  more  stable  charging  simulation. 

Nascap-2k  uses  a  high-order,  finite -element  representation  for  the  electrostatic  potential  that 
ensures  electric  fields  are  strictly  continuous  throughout  space.  The  electrostatic  potential  solver 
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uses  a  finite  element/conjugate  gradient  technique  to  solve  for  the  potentials  and  fields  on  the 
spacecraft  surface  and  throughout  the  surrounding  space.  Space  charge  density  models  presently 
include  Laplacian,  Linear,  Nonlinear,  Frozen  Ions,  Consistent  with  Ion  Density,  Full  PIC 
(Particle-in-Cell),  and  Hybrid  PIC. 

Particle  tracking  is  used  to  simulate  sheath  currents,  to  study  detector  response,  or  to  generate 
space  charge  evolution  for  dynamic  calculations.  Nascap-2k  generates  macroparticles  (each  of 
which  represents  a  collection  of  particles)  at  either  a  sheath  boundary,  the  problem  boundary,  at 
user- specified  locations,  or  throughout  all  space.  Particles  are  tracked  for  a  specified  amount  of 
time,  with  the  timestep  automatically  subdivided  at  each  step  of  each  particle  to  maintain 
accuracy.  The  current  to  each  surface  element  of  the  spacecraft  is  recorded  for  further  processing. 

Nascap-2k  User  s  Manual  is  the  primary  document  on  the  use  of  the  code.  This  document, 
Nascap-2k  Scientific  Documentation ,  describes  the  physics  and  numeric  models  used  in  the 
surface  charging,  potential  solution  and  particle  tracking  portions  of  the  code.  Section  2,  Physical 
Models,  describes  the  physics  models  included  in  the  code.  Section  3,  Numeric  Models  describes 
the  implementation. 

2  PHYSICAL  MODELS 

Nascap-2k  incorporates  physical  models  appropriate  to  both  tenuous  (e.g.,  GEO  orbit  or 
interplanetary  missions)  and  dense  (e.g.,  LEO  orbit)  plasma  environments.  The  code  solves  for 
environmentally-induced  time-dependent  potentials  on  spacecraft  surfaces,  the  electrostatic 
potential  about  the  object,  and  the  resulting  charged  particle  motion. 

Time-dependent  surface  potentials  are  computed  using  orbit  limited  currents  in  tenuous  plasma, 
tracked  currents  and/or  analytic  approximations  in  dense  plasma,  or  a  combination  in 
intermediate  environments. 

The  electrostatic  potential  solver  uses  a  conjugate  gradient  technique  with  flexible  boundary 
conditions  to  solve  Poisson’s  equation  for  the  potentials  and  fields  throughout  the  computational 
space.  Several  analytic  and  numeric  space  charge  density  models  are  available. 

Particle  tracking  is  used  to  study  sheath  currents,  to  study  detector  response,  to  generate  steady- 
state  charge  densities,  and  to  generate  space  charge  evolution  for  dynamic  calculations. 

2.1  Surface  Charging  from  Orbit  Limited  Currents 

Spacecraft  surface  charging  is  the  accumulation  of  charge  on  spacecraft  surfaces.  As  illustrated 
in  Figure  1,  several  different  current  components  can  contribute  to  the  charging.  [6, 7, 8,9]  First 
there  are  the  incident  charged  particles.  The  impact  of  high-energy  electrons  and  ions  causes  the 
ejection  of  secondary  electrons  into  space.  Incident  electrons  can  also  be  reflected  as  backscatter. 
In  sunlight  solar  ultraviolet  radiation  generates  low-energy  photoelectrons,  which  are  ejected 
from  the  surface.  Due  to  the  low  mass  of  electrons,  the  incident  electron  current  is  generally 
higher  than  the  ion  current.  By  itself,  this  would  lead  to  negative  charging.  However,  with  the 
inclusion  of  the  other  terms,  all  of  which  contribute  positive  current,  the  net  current  can  be  either 
positive  or  negative. 
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secondary,  backscattered, 
photo  electrons 


Figure  1.  High  Negative  Potentials  can  Result  from  the  Accumulation  of  Charge  on 

Spacecraft  Surfaces 

The  rest  of  the  spacecraft  influences  the  potential  and  local  electric  field  of  each  surface.  Each 
conductor  and  each  insulating  spacecraft  surface  collects  current  from  the  plasma  and  is 
capacitively  and  resistively  coupled  to  the  other  conductors  and  other  surfaces.  Additionally, 
some  spacecraft  have  current  sources,  such  as  electron  and  ion  guns,  that  contribute  to  the  net 
current.  Nascap-2k  uses  the  spacecraft  geometry,  surface  materials,  plasma  environment,  and  sun 
intensity  and  direction  to  calculate  the  time  history  of  the  surface  potentials,  electric  fields,  and 
fluxes. 

If  the  net  current  is  initially  negative  or  positive,  the  exposed  surface  begins  to  acquire  a  negative 
or  positive  potential  respectively.  As  the  magnitude  of  the  potential  increases,  the  net  current  is 
attenuated,  until  it  eventually  approaches  zero  and  the  surface  potential  reaches  its  equilibrium 
value.  Negative  potentials  as  high  as  ten  kilovolts  have  been  observed  during  geosynchronous 
substorms  and  over  one  kilovolt  during  auroral  precipitation  events.  Because  positive  currents 
generally  result  from  emission  of  low  energy  secondary  electrons  and  photoelectrons,  the 
positive  potentials  that  can  be  attained  are  relatively  modest.  As  soon  as  the  surface  reaches  a 
potential  greater  than  that  of  the  emitted  electrons  they  are  re-attracted  to  the  surface,  so  their 
contribution  to  the  current  is  suppressed.  The  equilibrium  potential  is  determined  by  the 
suppression  of  low  energy  emission  due  to  the  surface’s  own  electric  field.  A  similar  suppression 
effect  may  occur  due  to  the  electric  fields  of  neighboring  negative  charged  surfaces.  This  effect  is 
the  primary  reason  that  accurate  spacecraft  surface  charging  calculations  are  inherently  three 
dimensional. 

Positive  potentials  up  to  100  V  are  possible  in  environments  where  the  plasma  current  is  small 
compared  with  the  photo  current  (outer  magnetosphere  and  in  the  solar  wind).  In  these 
environments,  the  small  incident  plasma  current  is  balanced  by  escape  of  the  high  energy  tail  of 
the  photoelectron  energy  distribution,  even  though  the  bulk  of  the  emitted  photoelectrons  return 
to  the  surface. 

Nascap-2k  models  the  spacecraft  as  a  collection  of  surfaces  that  are  either  conductive  or  are  a 
dielectric  film  over  a  conductive  substrate.  The  conductive  surfaces  and  substrates  can  be  at  the 
same  potential  or  biased  with  respect  to  each  other. 
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2.1.1  Plasma  Environment 


In  the  orbit  limited  approximation  of  the  current  incident  on  a  surface,  the  flux  is  calculated  as 
that  incident  on  a  sphere  of  the  same  potential  and  depends  only  on  the  surface  potential,  and 
optionally  an  angle  with  respect  to  the  direction  of  flow.  This  approximation  applies  if  any 
charged  particle  far  from  the  spacecraft  can  reach  the  surface,  that  is,  there  are  no  excluded 
orbits.  For  spacecraft  this  is  generally  true  if  the  Debye  length  is  long  compared  with  the 
spacecraft  size  and  the  spacecraft  is  essentially  convex  with  no  self-shading.  “Surface  Charging” 
calculations  with  “Analytic  Currents”  in  “Geosynchronous,”  “Auroral,”  and  “Interplanetary” 
environments  use  the  orbit  limited  approximation,  and  “Surface  Charging”  calculations  with 
“Tracked  Ions  and  Analytic  Electron  Currents”  in  the  “Auroral”  environment  use  this 
approximation  for  electrons. 

The  analytic  current  density  to  a  surface  at  potential  (j>  is  given  by  integrating  over  the 
distribution  function. 


j  =  q 


~  2Ee2 

f  E  1 

nr 

J 

,E±<|>, 

f  (E  ±  (j))dE  =  q 


E  ±  (j) 


F(E  ±  (j))dE 


(1) 


Where  the  integration  variable  E  represents  the  energy  at  the  surface,  the  integration  limit  L  is  0 
for  the  repelled  species  and  l(j)l  for  the  attracted  species,  the  fraction  in  parentheses  represents  the 
orbit  limited  enhancement  or  reduction  of  current,  the  upper  sign  is  for  ions  and  the  lower  sign  is 
for  electrons,  f(E)  is  the  distribution  function  at  infinity,  and  F(E)  is  the  flux  at  infinity. 
Depending  on  the  model  selected  the  flux  is  given  by  one  of  the  following  formulas,  in  which  the 
particle  energy  E  and  the  temperature  9  are  measured  in  electron  volts,  and  the  surface  potential 
4>  is  measured  in  volts. 
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This  distribution  is  used  to  model  auroral  electrons. 
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where  H  is  the  Heaviside  step  function  and  El  and  Eu  are  upper  and  lower  energy  limits. 
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Convected  Maxwellian 


This  distribution  is  used  to  model  a  flowing  plasma. 
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(6) 


where  U  is  the  velocity  of  the  plasma  in  the  spacecraft  frame  of  reference  and  x  is  the  angle 
between  the  flow  vector  and  the  incident  velocity  at  infinity,  which  can  be  related  to  the  velocity 
at  which  the  particle  strikes  the  surface. 

2.1.2  Secondary  and  Backscattered  Current 

The  secondary  and  backscattered  currents  are  given  by 


Js.b  =  q 


Ys'b(E) 


E  ±  (j) 


F(E  ±  (j))dE 


(7) 


where  Ys  b  is  the  number  of  secondary  or  backscattered  (respectively)  electrons  per  incident 
charged  particle  (after  accounting  for  the  angular  distribution  of  the  incident  particles)  and  F  is 
the  flux  at  infinity  of  the  appropriate  incident  charged  particle. 


We  take  the  spectrum  of  emitted  secondary  electrons  to  be  a  Maxwellian  characterized  by  a 
temperature  of  2  eV.  It  follows  that  the  effective  secondary  electron  current  from  surfaces  at 
positive  potentials  is  reduced  by  a  factor  of  exp(-(j)/2)  relative  to  Equation  (7).  The  current  can 
also  be  reduced  by  an  electron-attracting  (E  •  n  >  0)  electric  field. 


2.1.3  Photocurrent 


For  surfaces  at  negative  potentials,  the  photocurrent  depends  on  the  surface  material,  angle  of 
incident  sunlight,  and  distance  from  the  sun.  As  with  secondaries,  the  photoelectron  spectrum  is 
by  default  taken  to  be  a  Maxwellian  characterized  by  a  temperature  of  2  eV,  so  that  current 
escaping  from  surfaces  at  positive  potentials  is  reduced  by  a  factor  of  exp(-cj)/2). 

Spacecraft  in  the  solar  wind  normally  charge  to  positive  potential  (a  few  volts  to  tens  of  volts)  to 
balance  the  emitted  photoelectron  current  with  a  very  low  incident  current  of  ambient  electrons. 
Unlike  geosynchronous  substorm  charging,  where  the  spacecraft  goes  negative,  determining  how 
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positive  the  spacecraft  charges  requires  knowledge  of  the  high-energy  portion  of  the 
photoelectron  spectrum.  When  the  user  provides  a  tabular  spectrum,  the  tabular  spectrum  is  used 
to  determine  the  escaping  photocurrent  rather  than  the  default  exp(-(j)/2). 

2.1.4  Other  Current  Sources 

Conduction  Current  through  Insulators 

For  the  case  of  interest,  in  which  the  deposition  thickness  is  short  compared  with  the  thickness  of 
the  insulator  layer,  the  electric  field  is  spatially  uniform  throughout  the  layer  and  current  is 
conducted  through  it.  (Note:  the  range  of  a  40  keV  electron  in  aluminum  is  -15  microns 
(~  0.0006  inches).)  For  insulating  surface  elements,  the  current  through  the  layer  to  the 
underlying  conductive  substrate  is  given  as  a  function  of  the  conductivity  and  thickness  of  the 
layer  by  Ohm’s  Law: 


_  Acj> 

J  conductor  . 

Pd  (8) 

where  Acj)  is  the  difference  in  potential  between  the  surface  and  the  conductive  substrate  and  p 
and  d  are  the  resistivity  and  thickness  of  the  layer. 

Surface  Conductivity 

Nascap-2k  accounts  for  current  flowing  along  surfaces.  Current  flows  between  insulating  surface 
elements  of  a  common  material  and  with  a  common  edge,  for  transport  over  a  wide  expanse  of 
such  material.  The  surface  element  to  surface  element  conductivity  is  proportional  to  the  length 
of  the  common  edge  and  inversely  proportional  to  the  sum  of  the  centroid-to-edges  distances. 

For  a  surface  element  i  having  a  common  edge  with  j,  R;  =  (d/l.j)  p,  where  p  is  the  surface 
resistivity,  di  is  the  distance  from  the  centroid  of  surface  element  i  to  the  center  of  the  common 
edge,  and  lij  is  the  length  of  the  common  edge.  The  surface  conductance  across  the  edge  between 
adjoining  surface  elements  i  and  j  is 


Conductivity  = - .  (9) 

Ri  +  Rj 

A  surface  can  be  grounded  by  a  strip  at  an  element  edge  (grounding  edge)  or  by  a  circular  contact 
located  at  a  node  (grounding  node). 

Charged  Particle  Emission 

Nascap-2k  also  allows  for  the  inclusion  of  the  current  from  charged  particle  emission  from  a 
conductor  (e.g.,  an  electron  gun).  This  capability  is  accessed  by  editing  the  script. 

2.1.5  Timescales 

Charged  particles  with  energies  below  50  keV  penetrate  less  than  a  hundred  microns  into  the 
spacecraft  skin.  While  the  time  for  the  entire  spacecraft  to  achieve  net  current  balance  is  typically 
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several  milliseconds,  the  time  for  each  surface  to  achieve  its  own  equilibrium  potential  is 
thousands  of  times  longer.  The  development  of  differences  between  the  potentials  of  different 
surfaces  is  referred  to  as  differential  charging. 

The  rate  of  absolute  charging  is  determined  by  the  capacitance  of  the  spacecraft  to  infinity,  while 
the  rate  of  differential  charging  is  determined  by  the  capacitance  across  the  insulator  layer.  The 
rate  of  change  of  the  potential  on  a  sphere  is  given  by 

d(b  1  I  JR 

—  = - =  —  (10) 

dt  47T£0  R  £0 

where  I  is  the  current  to  the  sphere,  J  is  the  current  density  to  the  sphere  surface,  and  R  is  the 
sphere  radius.  For  a  1 -meter  sphere  and  a  net  current  of  1  pA  m~  ,  the  sphere  begins  to  charge  at  a 
rate  of  100  kV  s’1.  As  the  sphere  charges  toward  its  steady-state  (zero  net  current)  potential,  the 
net  current  decreases  toward  zero  as  electrons  are  repelled.  Therefore,  the  charging  rate  decreases 
with  time. 


The  rate  of  change  of  the  difference  in  potential  across  a  dielectric  layer  is  given  by 

dcj)  _  Jd 
dt  £0 


(11) 


where  d  is  the  thickness  of  the  layer.  A  100  pm  layer  with  the  same  incident  current  density 
charges  at  a  rate  of  10  V  s’1.  In  this  representative  example,  differential  charging  takes  place  104 
times  more  slowly  than  absolute  charging. 

2.1.6  Circuit  Model 


Figure  2  shows  a  circuit  diagram  for  a  spacecraft  with  one  insulating  surface  element  and 
exposed  conducting  surfaces.  The  widely  differing  capacitances  of  surfaces  to  infinity,  CA  and  Cs 
and  of  the  surface  to  spacecraft  ground,  CAs,  make  this  a  complex  numeric  problem. 

CAS  =££0-«SxlO'6  Farad  (12) 

d 


C ,  «  C 


4  71  £q  R 


f  s  ) 

rsl 

v4;tR2  j 

X 10  “Farad 


(13) 


where  £,  d,  and  S  are  the  dielectric  constant,  thickness  (m),  and  surface  area  (m2)  of  the 
insulating  surface  element,  and  R  is  the  radius  of  a  sphere  approximating  the  size  of  the 
spacecraft.  The  potentials  as  a  function  of  time  are  computed  using  implicit  time  integration  of 
the  charging  equations,  which  relate  the  derivative  of  the  potential,  <D,  with  time  to  the  current,  I. 


cAoA+cAS(<i»A-<i»s)=iA 
-cAS(d»A  -<s>s)+cs<bA  =iA 
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Insulating  surface 


Figure  2.  Circuit  Model  of  a  Spacecraft  with  One  Insulating  Surface  Element 

The  multisurface  problem  is  solved  by  linearizing  the  currents  and  inverting  the  matrix. 

cd>  =  l(<I>)  (15) 


2.2  Material  Properties 

The  properties  of  the  spacecraft  surface  materials  determine  the  charging  rates  and  equilibrium 
surface  potentials.  Nascap-2k  uses  the  models  of  secondary  electrons  due  to  incident  electrons, 
secondary  electrons  due  to  incident  ions,  backscattered  electrons,  and  photoelectrons  originally 
developed  for  NASCAP/GEO.  [10] 

2.2.1  Secondary  Electron  Emission  Due  to  Electron  Impact 


Secondary  electrons  are  those  emitted  from  a  surface  with  energies  below  50  eV  in  response  to 
energy  deposited  near  the  surface  by  incident  electrons.  Their  energy  distribution  is  usually 
peaked  below  10  eV.  The  secondary  yield,  Yee,  is  the  ratio  of  primary  to  secondary  electron 
current. 


Y„ 


emitted  secondary  current  due  to  electron  impact 
primary  electron  current 


A  typical  curve  is  shown  in  Figure  3. 


(16) 
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Figure  3.  Electron  Secondary  Yield  as  a  Function  of  Incident  Energy 

The  secondary  electron  emission  yield  can  be  calculated  using  the  empirical  formula  [11]: 


Y«(V)  =  C[ 

J( 


rR 

c 

dE 

Jo 

dx 

—ax  cos 


^dx 


(17) 


where  x  is  the  path  length  of  penetration  of  a  primary  electron  beam  into  the  material,  R  is  the 
“range,”  or  maximum  penetration  length,  and  \\t  is  the  angle  of  incidence  of  the  primary  electron. 

This  equation  is  based  on  a  simple  physical  model  [12]: 

1.  The  number  of  secondary  electrons  produced  by  the  primary  beam  at  a  distance  x  is 
proportional  to  the  energy  loss  of  the  beam  or  “stopping  power”  of  the  material, 
dE 
dx  ' 


2.  The  fraction  of  the  secondaries  that  migrate  to  the  surface  and  escape  decreases 

exponentially  with  depth  (f  =  c~axcosvl/  )  Thus  only  secondaries  produced  within  a  few 
multiples  of  the  distance  1/a  (the  depth  of  escape)  from  the  surface  contribute  significantly 
to  the  observed  yield. 

The  stopping  power  for  incident  electrons  of  energy  E  is  related  to  the  range  of  these  electrons 
through  the  equation 


S(E) 


dE 

dR 

dx 

dE 

(18) 


The  usual  formulation  for  the  range  is  that  it  increases  with  the  energy,  E  of  the  incident 
electrons  in  a  way  that  approximates  a  simple  power  law  [13]: 
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R=bEq 


(19) 


where  q  is  a  bit  less  than  2  for  electrons  with  energy  <  300  keV.  Because  the  primary  beam  loses 
energy  as  it  passes  through  the  material,  E  and  S(E0,  x)  depend  on  the  path  length  x,  where  E0  is 
the  initial  electron  energy.  The  stopping  power  can  be  written  as 


S(E0,x) 


dE 

dR 

-X_  1 

f  b  ^ 

dx 

dE 

qb 

(R-xJ 

(20) 


Figure  4  shows  S(E0,  x)  plotted  against  x  for  several  values  of  E0.  Inspection  of  Figure  4  and  the 
equation  for  S(x)  illustrates  the  following  points: 

1.  S(E0,  x)  increases  with  x,  slowly  at  first,  before  reaching  a  singularity  as  x  approaches  R. 

2.  The  initial  value  of  S(E0,  x)  decreases  with  increasing  initial  energy  E0. 

Both  of  these  observations  are  due  to  the  decrease  in  electron-atom  collision  cross-section  with 
increasing  energy. 


Path  Length,  x 


"  1/a  *■ 

Depth  of  Escape 


(a) 


(b) 


Figure  4.  Energy  Deposition  Profiles  of  Normally  Incident  Primary  Electrons  for  Four 
Incident  Energies  E° ,  E^ ,  E° ,  and  E”  and  Corresponding  Yield  Curve 

The  yield  is  only  sensitive  to  the  details  of  the  stopping-power  depth-dependence  for  initial 
energies  with  ranges  of  the  same  order  as  the  escape  depth,  R  ~  1/a  (i.e.,  about  the  maximum  of 
the  yield  curve).  For  lower  energies,  R  «  1/a,  essentially  all  the  primary  energy  is  available  for 
detectable  secondary  production,  leading  to  a  linear  increase  in  yield  with  increasing  Ec.  At 
higher  energies,  where  R  »  1/a,  S(E0,x)  remains  almost  constant  over  the  depth  of  escape. 
Therefore,  along  with  S(E0,x),  the  yield  decreases  as  E0  increases. 
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Taking  this  into  account,  the  stopping  power  can  be  approximated  by  a  linear  expansion  in  x, 
about  x  =  0. 


dF, 

f  dR  ^ 

-l 

^d2R  1 

f  dRl 

+ 

dx  ”4 

IdEj 

(21) 


The  simple  power  law  does  not  adequately  describe  the  available  experimental  data.  The  range 
has  different  exponents  for  low  energy  and  high  energy.  A  biexponential  range  law  with  four 
parameters  bi,  b2,  qi,  qofits  the  experimental  data  better. 


R=b,E’i +b2E’2  (22) 

Figure  5  shows  CSDA  range  values  for  Aluminum  for  energies  from  20  eV  to  100  keV.  Note  the 
drastic  change  in  behavior  near  300  eV.  The  biexponential  formula 
R  =  142.63E0'2335  +  245. 15E1 7269  fits  well  throughout  the  range  plotted. 


Electron  Energy  (keV) 


Figure  5.  CSDA  Range  for  Electrons  in  Aluminum.  Low  Energy  Values  (Blue  Diamonds) 
From  Ashley  etal.  (Reference  14).  High  Energy  Values  (Red  Squares)  From  NIST  ESTAR 

Database.  Dashed  Line  Is  Biexponential  Fit. 

For  materials  where  no  suitable  data  is  available,  a  monoexponential  form  can  be  generated 
using  Feldman’s  empirical  relationships  [13],  connecting  b  and  n  to  atomic  data. 


b  =  250w/pmaasZ‘>/2 

(23) 

q  =  1 .2/(1  -0,29 log |o  Z) 

(24) 

National  Institute  of  Standards  and  Technology  Electron  STopping  power  And  Range  database  available  at 
http://www.nist.gov/pml/data/star/index.cfm 
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where  W  is  the  atomic  or  molecular  weight  of  the  material,  Z  is  the  atomic  number,  and  pmass  is 
the  mass  density  in  gm  cm'3.  Then  E  is  in  keV  and  R  is  in  angstroms.  The  stopping  power  is  then 
obtained  indirectly  with  the  above  equation.  Theoretical  estimates  of  the  stopping  power  for  a 
number  of  materials  are  available  from  Ashley  et  al.  [14]  Comparison  of  these  values  with  those 
implied  by  the  range  data  showed  significant  discrepancies,  particularly  for  those  materials  fit 
using  Feldman’s  fonnula.  The  best  approach  is  to  fit  the  four  parameters  in  the  equation  for  R 
directly  to  the  stopping  power  data. 

S  =  (q1b1Ech-1+q2b2Ecl2-1)“1  (25) 

Burke  et  al.  [15]  propose  a  relationship  between  secondary  emission  due  to  electrons  below  1 
keV  and  secondary  emission  due  to  higher  energy  particles,  including  gamma  radiation  from 
Co60. 


As  is  evident  from  equation  (16),  the  secondary  electron  yield  depends  on  angular  distribution  of 
the  incident  electrons.  For  low  energy  electrons,  all  the  energy  is  deposited  near  the  surface 
regardless  of  incident  angle,  so  the  yield  is  independent  of  incident  angle.  For  high  energy 
electrons  the  stopping  power  is  constant  in  the  near-surface  region,  so  the  yield  is  proportional  to 
the  path  length  near  the  surface,  or  inversely  proportional  to  cos\|/.  It  follows  that  the  yield  for 
high  energy  isotropic  electrons  is  double  that  for  normally  incident  electrons  of  the  same  energy. 
The  isotropic  yield  curve  makes  a  transition  from  being  equal  to  the  normal  incidence  curve  at 
low  energy  to  double  at  high  energy.  The  maximum  yield  for  isotropic  electrons  is  somewhat 
higher  than  for  normal  electrons,  and  occurs  at  somewhat  higher  energy. 

2.2.2  Secondary  Electron  Emission  Due  to  Ion  Impact 

Secondary  emission  of  electrons  due  to  ion  impact  can  be  treated  in  a  way  similar  to  that  for 
electron  impact.  The  yield  is  given  by 


(26) 


The  stopping  power  is  assumed  to  be  independent  of  path  length  x  over  the  thickness,  d,  of  the 
sample.  At  low  energies  the  stopping  power  is  proportional  to  the  velocity,  and  at  high  energies  it 
is  inversely  proportional  to  the  velocity  [16]. 


dE  _  [3E1/2 

dx  1+ E/E, 


(27) 


Emax  is  the  energy  at  the  maximum  in  the  yield  curve.  This  is  approximately  50  keV  for  most 
materials. 


As  the  CSDA  range  is  at  least  comparable  to  the  mean  depth  of  secondary  electron  production 
for  protons  above  1  keV  (140  angstroms  of  aluminum),  we  consider  the  proton  yield  to  vary 
inversely  with  the  cosine  of  the  incident  angle,  so  that  the  isotropic  incidence  yield  is  double  the 
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normal  incidence  yield.  This  is  usually  an  adequate  treatment,  because  the  ion  fluxes  are  almost 
always  low,  and  the  effect  of  secondary  emission  is  to  mildly  augment  the  ion  flux,  rather  than 
canceling  it  as  is  the  case  with  incident  electrons. 

The  secondary  emission  properties  due  to  the  impact  of  ions  other  than  protons  are  assumed  to  be 
identical  to  the  proton  values  for  the  same  energy. 


The  secondary  electron  yield  curve  for  protons  incident  on  aluminum  is  shown  in  Figure  6.  As 
can  be  seen  by  comparison  with  Figure  7,  the  yield  closely  follows  the  stopping  power. 
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Figure  6.  Secondary  Electron  Emission  by  Aluminum  for  Proton  Impact  at  Normal 
Incidence;  Experimental  Points  as  Indicated  in  the  References  [17, 18, 19,  20] 
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Figure  7.  Stopping  Power  for  Protons  on  Aluminum  (from  NIST  PSTAR  Database1) 
2.2.3  Backscatter 

Backscattered  electrons  are  those  emitted  from  the  surface  with  energies  above  50  eV.  Their 
energy  distribution  is  usually  peaked  close  to  the  primary  incident  energy  and  they  may  be 
considered  as  reflected  electrons. 


Nascap-2k  uses  a  backscattering  theory  [21]  based  on  that  of  Everhart  [22]  as  extended  by 
McAfee.  [23]  It  assumes  a  single  scattering  in  accordance  with  the  Rutherford  cross-section  and 
the  Thomson- Widdington  slowing  down  law, 


dE 

dx 


oc  E 


(28) 


(valid  for  most  metals  for  E  >  10  keV).  For  normal  incidence  the  backscattering  coefficient  is 
given  by 


(29) 

where  a  is  taken  to  be  0.0375  Z  and  where  Z  is  the  atomic  number  of  the  material.  This 
expression  matches  the  experimental  data. 

The  large-angle  scattering  theory,  together  with  Monte  Carlo  data  and  experiments  by  Darlington 
and  Cosslett,  [24]  indicate  that  the  angular  dependence  of  backscattering  is  well-described  by 

Tl(v)  =  Bo  exp(Tll  (l  — cos\|/))  (30) 


1  National  Institute  of  Standards  and  Technology  Proton  Stopping  power  And  Range  database  available  at 
http://www.nist.gov/pml/data/star/index.cfm 
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where  the  value  of  rp  is,  within  the  uncertainty  in  the  data,  what  would  be  obtained  by  assuming 
total  backscattering  at  glancing  incidence,  r|L  =  -log  r|  .  The  net  albedo  for  an  isotropic  flux  is 
then 


A  _2l-TTo(l-lQgT|o) 

(iogTlo)2 


(31) 


As  the  energy  is  decreased  below  10  keV,  the  backscattering  increases.  Data  cited  by  Shimizu 
[25]  indicate  an  increase  of  about  0.1,  almost  independent  of  Z.  This  component  of 
backscattering  can  be  approximated  by 


8r|0  =0.1exp[-E/5  keV] 

At  very  low  energies,  the  backscattering  coefficient  becomes  very  small  and,  below  50  eV, 
becomes  zero  by  definition,  as  all  emitted  electrons  are  considered  secondaries.  This  can  be 
taken  into  account  by  a  factor  of 


H(E-50eV) 


log  20 


log 


E  ^ 

50eV  j 


(33) 


where  H  is  the  Heaviside  step  function.  The  formula  for  energy-dependent  backscattering, 
incorporating  these  assumptions,  is  then 


bo  = 


H(l-E)H(E-0.05) 


'  1  A 

log  20 


log 


0.05 


+  H(E-l) 


j-E/5 

"To" 
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^2A0.037Z^ 


(34) 


where  energies  are  measured  in  keV. 

2.2.4  Photoemission 


Photoelectrons  are  those  ejected  from  the  surface  due  to  the  solar  ultraviolet  radiation.  Usually 
the  quantity  known  is  the  yield,  or  number  of  electrons  emitted  for  a  surface  normally  exposed  to 
the  solar  spectrum,  at  an  “earth  distance”  from  the  sun.  The  photocurrent  from  a  surface  exposed 
to  the  sun  at  an  angle  v|/sun  is  given  by  the  formula 

i  photo  =  (Area  exposed)Ysun  cos  y  sun  (35) 


This  assumes  that  the  yield  per  photon  is,  on  average,  independent  of  \|/sun. 

A  surface  element  is  taken  to  be  sunlit  if  its  centroid  is  sunlit,  i.e.,  if  a  vector  from  the  centroid  in 
the  sun  direction  does  not  intersect  any  other  surface  element. 
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2.3  Surface  Charging  from  Tracked  Currents  and/or  Current  Balance 

In  short  Debye  length  environments,  such  as  encountered  in  low-Earth  orbit,  the  orbit-limited 
approximation  greatly  overestimates  currents  of  attracted  particles.  This  occurs  because 
potentials  that  decay  faster  than  inverse  square  (sheath  potentials)  lead  to  an  effective  potential 
barrier  to  attracted  charged  particles  when  angular  momentum  is  taken  into  account.  Nascap-2k 
can  compute  the  charging  of  spacecraft  surfaces  using  tracked  currents,  either  by  themselves  or 
in  conjunction  with  analytically  computed  currents.  Calculations  of  typical  interest  include: 

1.  Floating  potential  and/or  differential  charging  of  a  negative  spacecraft,  whose  ion  current 
may  be  enhanced  and  asymmetrized  due  to  ram/wake  effects. 

2.  Floating  potential  of  a  biased  spacecraft  that  has  both  ion  and  electron  sheaths.  The 
electron  current  may  be  additionally  limited  by  magnetic  field  effects. 

3.  Dynamic  potential  of  a  spacecraft  with  time-dependent  bias. 

In  “Surface  Charging”  calculations  with  “Tracked  Ion  and  Analytic  Electron  Currents” 
(“Auroral”  environment  only),  the  analytic  electrons  are  high  energy  electrons  computed  from 
the  Fontheim  formula.  In  “Time  Dependent  Plasma”  calculations,  the  optional  electron  thermal 
current  is  given  by 


Ajth  exp(c|)/e) 


Aith 


l4 


if  4>  <  0 
if  c()  >  0 


(36) 


where  jth  =  ne^^— 

2.4  Potentials  in  Space 

Nascap-2k  uses  a  finite  element  method  to  solve  Poisson’s  equation  -s0V2c|)  =  p(x)  for 
electrostatic  potentials  throughout  space.  Nascap-2k  uses  the  variational  form 


o=A{ 

dV 

^  IJ 

L2  £o  J 

+  jc|)Vc|)»dS 


(37) 


The  first  term  in  the  volume  integrand  corresponds  to  the  Faplacian  operator;  the  second  term  is 
the  space  charge  contribution.  As  surface  potentials  are  held  fixed  when  solving  for  potentials  in 
space,  the  surface  term  appearing  in  equation  (36)  is  not  used  in  the  potential  solution,  but  may 
be  used  later  to  determine  the  normal  component  of  electric  field  (and  thus  the  surface  charge) 
for  each  surface  element. 

The  solution  is  subject  to  fixed  potential  boundary  condition  on  the  spacecraft  surfaces  and  either 
zero  potential  or  monopole  potentials  on  the  grid  boundary. 

A  number  of  models  of  the  charge  density  are  available  to  study  different  phenomena. 
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In  the  following,  the  symbol  A,d  is  the  plasma  Debye  length  and  g  is  the  plasma  density  reduction 
factor  computed  by  geometric  wake  model  (0<g<l). 

Laplace.  The  Laplacian  space  charge  option  specifies  that  the  charge  density  is  zero  and  is 
appropriate  for  very  low  density  plasmas. 


8 


O 


(38) 


i.e.,  charge  exists  only  on  object  surfaces  and  external  boundaries,  as  determined  by  the 
boundary  conditions.  “Space  charge”  iterations  may  still  be  required,  however,  due  to  the 
treatment  of  surface  electric  fields. 


Linear  (Debye  Shielding).  The  Linear  space  charge  option  solves  the  Helmholtz  or  Debye  - 
Huckel  equation.  This  model  is  very  rarely  used. 


-vV  = 


(39) 


Nonlinear.  The  Nonlinear  space  charge  model  is  appropriate  for  most  low-Earth-orbit  type 
plasmas.  It  accounts  for  space  charge  acceleration  and  convergence  in  a  manner  based  on 
spherical  collection  (Langmuir-Blodgett  [26]  problem).  Poisson’s  equation  is  solved  with  space 
charge  given  by: 


p  _ 

(  max(l,  c(cf>. 
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(40) 


This  model  smoothly  interpolates  between  Debye  screening  at  low  potential  and  an  accelerated 
distribution  with  particle  convergence  at  high  potentials.  The  convergence  factor  is  computed  in 
terms  of  local  information  and  problem  parameters.  The  convergence  model  was  developed  by 
numerically  solving  the  Langmuir-Blodgett  problem  for  collection  by  a  high-voltage  sphere  and 
fitting  the  result  to  an  analytic  form.  An  excellent  fit  was  found  with 


(R,h/r)I=2.29||E|XD/e|‘“2|0/*rJ”’ 

(41) 

C(f  |E|)=  min((R  „  /r)2 ,3.545|<|>/e| K ) 

(42) 

where  E  is  the  local  electric  field. 

Frozen  Ion.  The  frozen  ion  formulation  is  intended  for  short  timescale  (typically  sub¬ 
microsecond)  problems  for  which  it  is  a  good  approximation  to  assume  that  ions  remain 
stationary  and  at  ambient  density  (“ion  matrix”  approximation),  but  electrons  achieve  barometric 
equilibrium.  This  approximation  is  intended  primarily  for  negative  potentials,  for  which  the 
charge  density  is  zero  at  zero  potential  and  at  negative  potential  (electrons  repelled)  the  charge 
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density  becomes  positive  and  asymptotically  reaches  the  ion  density.  For  computational 
convenience  we  extend  the  charge  density  linearly  to  positive  potentials,  as  an  exponential 
increase  in  electron  charge  would  be  unphysical. 

p  =  ne(l-e*/e)  §<0  (43) 

P =  -ne  q  (t)>0  (44) 

Considerations  for  implementing  this  formulation  stably  on  a  coarse  mesh  are  discussed  in 
Section  3.4.5. 

Full  Trajectory  Ions.  Ion  densities  are  calculated  from  steady-state  ion  trajectories.  Electrons 
are  barometric.  This  approximation  is  used  to  model  charge  density  for  a  moving  spacecraft 
when  current  collection  in  the  wake  is  important  and  the  neutral  wake  approximation  is 
inadequate — such  as  when  high  potential  surfaces  are  in  the  wake. 


~  =  ~  (l -  exp  (((j)  -  (j)b  )/9)) 

eo  £0 

=01n( 


^Pi A 


yen; 


(45) 


Here,  p,  is  calculated  by  tracking  ions. 

Plume  Ion  Density.  Ion  densities  are  initially  computed  from  an  imported  plume  map  file. 
Electrons  are  barometric.  The  quantity  p/e0  is  computed  using  the  same  formula  as  for  full 

trajectory  ions,  where  pi  is  computed  by  summing  the  contributions  from  the  thruster  plumes 
instead  of  computed  by  tracking  macroparticles.  After  the  initial  iteration,  pi  is  computed  by 
summing  the  main  beam  contribution  from  the  thruster  plumes  and  the  charge  exchange 
contribution  from  generating  and  tracking  charge  exchange  ions. 

Hybrid  PIC.  This  algorithm  is  used  for  timescales  (typically  sub-millisecond)  on  which  it  is 
practical  to  treat  ion  motion,  but  electrons  may  be  considered  in  barometric  equilibrium.  The 
total  charge  density  is  the  sum  of  the  ion  and  electron  charge  densities. 


P  =  Pe  +Pi 


(46) 


The  ion  density  is  computed  from  ion  macroparticles  that  move  in  the  local  electric  fields  at  each 
timestep.  The  electron  charge  density  is  given  by 


Pe 
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(47) 
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Full  PIC.  For  this  option,  it  is  assumed  that  both  the  electron  and  ion-charge  densities  were 
stored  during  particle  tracking. 


_  _  tracked 

P  =  P 


(48) 


2.5  Macroparticle  Creation  and  Tracking 

Particle  tracking  is  used  to  study  sheath  currents,  to  study  detector  response,  to  generate  steady- 
state  charge  densities,  and  to  generate  space  charge  evolution  for  dynamic  calculations.  The 
investigation  of  detector  response  is  applicable  in  both  tenuous  and  dense  plasmas.  However,  the 
other  phenomena  are  only  applicable  in  plasmas  with  Debye  lengths  no  more  than  the  spacecraft 
size. 


In  a  dense  plasma  with  high  potentials  the  disturbed  region  is  referred  to  as  a  sheath.  The  “sheath 
surface”  represents  a  sharp  demarcation  between  a  low  potential  exterior  region  containing 
neutral  “undisturbed  plasma,”  and  a  high  potential  interior  region  from  which  one  species  (ions 
for  positive  potentials;  electrons  for  negative  potentials)  is  excluded. 

In  a  dense  plasma,  the  total  current  and  its  distribution  over  the  spacecraft  surface  depend  on  the 
plasma  environment,  the  surface  potentials,  the  sheath  structure,  the  spacecraft  velocity,  and 
ambient  and  spacecraft  generated  magnetic  fields.  For  a  large  object  (many  Debye  lengths  in 
size)  at  low  potential  (at  most,  a  few  times  the  plasma  temperature),  Debye  screening  minimizes 
the  extent  of  the  object’s  influence  on  the  plasma.  Because  the  sheath  is  “thin,”  the  sheath  surface 
lies  near  the  spacecraft  surface.  The  current  attracted  to  each  surface  element  is  then  given  by  its 
area  times  the  plasma  thermal  current  of  the  appropriate  species.  For  an  object  (not  necessarily 
large)  at  high  potential,  however,  the  extent  of  perturbed  plasma  is  far  greater.  A  first  estimate  of 
the  sheath  thickness  is  given  by  the  Child-Langumir  [27,  28]  distance. 
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The  plasma  thermal  current  then  flows  through  an  area  much  larger  than  the  spacecraft  surface. 
The  current  may  be  focused  by  the  fields  within  the  sheath,  which  affects  the  distribution  of 
current  over  the  spacecraft  surface.  This  focusing  causes  the  increase  in  space  charge  density  due 
to  particle  acceleration  modeled  in  the  “Non-Linear”  space  charge  density  model. 

Macroparticles  are  generated  in  the  Particle  Generation  step  of  the  script  and  tracked  in  the 
Track  step  of  the  script.  Macroparticles  can  be  generated  at  a  well  defined  sheath  surface,  at  the 
boundary  of  the  computational  space,  throughout  the  volume  of  space,  or  at  user  specified 
locations.  As  macroparticles  are  tracked,  their  charge  density  in  the  volume  elements  they  pass 
through  may  or  may  not  be  retained,  depending  on  user  choices.  The  current  to  any  surface 
element  hit  by  the  macroparticles  is  retained. 

2.5.1  Generation  of  Macroparticles  at  Sheath  Surface 

The  attracted  species  diffuses  across  the  sheath  surface  at  a  rate  given  by  the  plasma  thermal 
current. 
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The  current  density  to  the  sheath  edge  is  computed  by  assuming  the  sheath  to  be  a  perfectly 
absorbing  spherical  surface. 

The  charged  particles  constituting  this  current  undergo  motion  consistent  with  the  electric  and 
magnetic  fields  inside  the  sheath.  The  sheath  surface  is  taken  as  the  equipotential  surface  at 
(j)  =  ±01n2.  This  choice  is  made  because  the  attracted  species  is  absorbed  by  the  sheath,  so  we 
have  only  the  inward  moving  component,  comprising  half  the  ambient  density.  The  repelled 
species,  whose  density  satisfies  n  =  na  exp(- 14>/ 0|) ,  must  also  be  at  half  the  ambient  density, 
leading  immediately  to  the  sheath  potential. 

Macroparticles  representing  this  sheath  current  can  be  generated  and  tracked  to  determine  the 
current  to  the  spacecraft  and  its  distribution.  The  initial  velocity  is  the  average  velocity  of  a 
charged  particle  crossing  the  sheath  boundary.  In  the  absence  of  spacecraft  motion  or  a  magnetic 

(2e0 

field,  this  is  v  =  -  in  the  direction  of  the  local  electric  field.  More  precisely,  for  a  sheath 

V  rtm 

segment  on  the  ram  side 
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where  the  plus  sign  is  for  ions  and  the  minus  for  electrons  and  vx  is  the  component  of  the 
spacecraft  velocity  along  the  local  electric  field.  In  the  absence  of  spacecraft  motion  and  for  a 
sheath  segment  in  the  wake, 


With  a  magnetic  field 
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(53) 


where  the  drift  velocity  term  vdrift 
square  root  negative. 


ExB 

B2 


is  not  included  if  it  would  make  the  argument  of  the 


This  model  is  appropriate  for  high  potentials.  When  either  the  thermal  distribution  of  particle 
velocities  or  the  spacecraft  velocity  is  important,  macroparticles  should  be  tracked  from  the 
boundary  of  the  computational  space. 
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2.5.2  Generation  of  a  Thermal  Distribution  at  the  Computational  Space  Boundary 

There  are  two  models  for  creation  of  macroparticles  at  the  boundary  of  the  computational  grid.  In 
both  models,  the  plasma  is  assumed  to  be  Maxwellian.  The  two  models  are  inconsistent  with 
each  other  and  should  not  both  be  used  together. 


Boundary  Injection 

In  the  simplest  model,  macroparticles  are  created  with  charge  equal  to  the  plasma  thermal  current 

(2e0 

times  the  area  times  the  time  interval,  q  =  jthAAT,  and  velocity  equal  to  , -  normal  to  the 


rcm 


boundary,  so  that  they  represent  a  density  of  n/2.  These  can  each  be  split  into  eight  outwardly 
moving  macroparticles,  to  approximate  a  thermal  distribution.  The  split  macroparticles  have  an 

feO 

additional  velocity  of  ±0.707  J —  along  each  axis  of  a  randomly  generated  coordinate  system. 


m 


The  splitting  is  done  in  the  plasma  frame  of  reference  in  order  to  simulate  the  correct  momentum 
and  energy  distribution  for  a  drifting  Maxwellian  when  transformed  back  to  the  spacecraft 
reference  frame. 


Boundary 

In  the  other  model,  macroparticles  that  represent  a  specified  fraction  of  the  thermal  distribution 
in  each  of  the  three  spatial  directions  are  created.  At  each  boundary  emission  point,  several 
macroparticles  are  created  that  sample  the  velocity  distribution  function  in  each  of  the  three 
spatial  directions.  The  user  specifies  which  fraction  of  the  distribution  function  each 
macroparticle  is  to  represent  and  the  code  computes  the  range  of  velocities  (v*’2y,z)  which  would 
give  the  specified  fractions,  tp. 
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Several  macroparticles  are  created,  each  with  velocity  the  sum  of  the  average  velocity  of  its 
fraction  of  the  distribution  and  the  ram  velocity  (negative  of  the  spacecraft  velocity). 
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2.5.3  Generation  of  Macroparticles  Throughout  the  Computational  Space 
Initialization 

Macroparticles  can  be  created  throughout  the  computational  space  to  represent  a  uniform 
distribution  of  charge.  The  weight  of  these  macroparticles  is  decreased  by  any  geometric  wake 
factor.  Each  macroparticle  can  be  split  into  eight  outwardly  moving  macroparticles  so  as  to 

approximate  a  thermal  distribution.  The  velocities  have  components  of  ±0.707  yeO/'rn  along 

each  axis  in  a  randomly  oriented  coordinate  system.  The  splitting  is  done  in  the  plasma  frame  of 
reference  in  order  to  simulate  the  correct  momentum  and  energy  distribution  of  a  drifting 
Maxwellian  when  transformed  back  to  the  spacecraft  reference  frame. 

Charge  Exchange 

Macroparticles  can  also  be  created  throughout  the  computational  space  to  represent  charge 
exchange  current.  Charge  exchange  ions  are  those  created  by  a  charge  transfer  collision  between 
a  fast  ion  and  a  slow  neutral  atom,  such  as  occurs  in  ion  thruster  plumes.  The  current  per  unit 
volume  associated  with  each  macroparticle  is  the  product  of  the  charge  exchange  cross  section, 
the  neutral  density,  and  the  main  beam  ion  current  density. 

J  cex  C ^ c ex  FI  n e li t ra b  bea m  V beam-  (56 


2e0  3448 

Each  component  of  the  initial  velocity  is  given  by  C  - : - ,  where  0.3448  is  4000  K  in 


m 


electron  volts  and  C,  is  a  random  number  between  -1  and  1. 

The  neutral  density  is  the  sum  of  the  un-ionized  propellant  from  the  thrusters,  the  gas  flowing 
through  the  neutralizers,  and  the  background  gas. 
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The  contribution  from  each  thruster  is  given  by  the  solid  angle  subtended  by  the  thruster  grid,  the 
flow  rate  and  the  temperature. 


unionized 
1  thrusters 


Flow 


neutral 

thruster 


71111  thruster 

V  thruster 


(58) 


The  contribution  from  each  neutralizer  is  given  by  the  flow  rate,  the  temperature,  the  angle 
between  the  neutralizer  axis  and  the  line  of  sight  from  the  neutralizer,  and  the  distance  from  the 
neutralizer. 
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2.5.4  Particle  Tracking 

Macroparticles  are  tracked  until  they  leave  the  computational  grid,  hit  a  surface,  exceed  the 
number  of  allowed  steps,  or  exceed  the  computational  tracking  time.  They  are  tracked  in  the 
local  electric  field  and  the  uniform  magnetic  field  using  a  third  order  energy  conserving 
algorithm. 

As  the  macroparticles  pass  through  volume  elements,  their  charge  and/or  current  can  be 
optionally  saved  to  volume  elements  and  nodes.  This  charge  and/or  current  can  subsequently  be 
used  in  the  space  charge  formulas  to  compute  potentials  in  space. 

2.6  Geometric  Wake 

The  exclusion  of  plasma  in  the  wake  of  a  rapidly  moving  spacecraft  is  also  important.  That 
portion  of  the  sheath  surface  that  lies  in  the  wake  region  has  a  reduced  current  passing  through  it. 
The  volume  of  space  in  the  wake  has  a  reduced  charge  density.  Conversely,  on  the  ram  side  of 
the  spacecraft,  there  is  enhanced  current  and  charge  density. 

A  spacecraft  creates  a  wake  when  its  velocity  is  comparable  to,  or  larger  than,  the  thermal 
velocity  of  the  ambient  ions.  In  low  Earth  orbit,  the  spacecraft  velocity  is  7500  ms"1.  The  thermal 
velocity  of  a  0.3  eV  oxygen  ion  is  only  about  2000  ms"1.  Thus  the  spacecraft  travels  a  few  of  its 
own  radii  before  ions  can  fill  in  behind  it.  The  electrons  fill  in  more  rapidly.  However,  the 
density  to  which  electrons  can  accumulate  is  limited  by  the  space  charge  of  the  electrons  already 
in  the  wake.  Thus,  except  in  regions  where  the  density  is  extremely  low,  the  electron  density  is 
only  very  slightly  greater  than  the  ion  density. 

For  each  point  in  space,  the  wake  module  calculates  what  fraction  of  the  thermal  distribution  of 
velocities,  (in  the  spacecraft  frame,  a  Maxwellian  displaced  by  the  spacecraft  velocity  vector)  is 
not  blocked  by  the  spacecraft. 

The  neutral  wake  density  is  given  by 


g(x,Q 


(60) 


where  fio(v)  is  the  unperturbed  distribution  function  for  a  drifting  Maxwellian,  and  g(x,Q)  has 
value  “0”  if  a  ray  starting  from  x  and  going  in  the  direction  Q  would  strike  the  spacecraft  and  “1” 
if  it  would  not. 

2.7  Plumes 

The  ion  density  and  velocity  due  to  an  ion  thruster  main  beam  can  be  read  into  Nascap-2k  and 
displayed.  Multiple  thrusters  (each  with  its  own  location  and  axis  direction)  can  be  specified, 
although  at  present  all  thrusters  have  the  same  type  of  plume.  Neutralizers  can  be  specified  as 
well.  The  plume  currents  are  used  in  the  computation  of  the  potentials  in  space  (“Plume  Ion 
Density”  space  charge  model)  and  charge  exchange  currents. 
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2.8  Transverse  Surface  Currents 


The  transverse  surface  current  algorithm  in  Nascap-2k  is  used  to  calculate  the  current  flowing 
along  an  antenna  element  or  other  object  in  response  to  time-varying  applied  potentials  and/or 
active  or  passive  current  sources.  The  model  accounts  correctly  for  both  local  capacitance  and 
incident  plasma  current.  The  change  in  charge  on  a  surface  element  during  a  timestep  is  that 
which  is  needed  to  accomplish  the  change  in  surface  electric  field,  which  is  obtained  from  the 
Nascap-2k  potential  calculation.  The  current  to  each  surface  element  consists  of  the  transverse 
surface  current  from  neighboring  surface  elements  and  the  current  provided  by  the  plasma,  which 
is  provided  by  the  Nascap-2k  particle-in-cell  (PIC)  simulation.  The  current  continuity  equation  is 
solved  using  a  pseudopotential  approach.  As  a  boundary  condition,  one  surface  element  of  each 
conductive  element  (i.e.,  surfaces  associated  with  a  given  Nascap-2k  conductor)  is  specified  as 
connected  to  the  biasing  power  supply,  and  the  solution  provides  the  required  current. 

The  basic  equation  to  be  solved  is  V  •  J  +  ^  =  0,  where  J  is  the  surface  current  density,  and  p  is 

the  surface  charge  density.  Note  that  the  solution  for  the  current  is  non-unique  to  the  addition  of 
any  divergence-free  current  field.  We  assume  that,  provided  appropriate  boundary  conditions  are 
implemented,  such  circulating  currents  are  not  of  concern  to  the  problem  being  solved. 

3  NUMERIC  MODELS 


This  section  describes  the  numeric  techniques  Nascap-2k  uses  to  implement  the  physics  models 

described  in  the  previous  section. 

•  The  orbit  limited  currents  are  computed  by  integrating  over  the  distribution  function.  An 
implicitized  version  of  the  Boundary  Element  Method  is  used  to  compute  the  evolution  of 
surface  potentials. 

•  The  finite  element  approach  is  used  to  solve  Poisson’s  equation  for  the  spatial  potential. 
Implicitization  and  charge  stabilization  are  used  for  stability.  The  finite  element  grid  uses 
strictly  continuous  electric  field  interpolants  for  accurate  particle  tracking. 

•  Macroparticle  tracking  is  used  for  computing  surface  and  volume  currents  and  volume 
densities.  Macroparticle  splitting  is  used  for  improved  accuracy. 

3.1  Environment  Integrals 

In  section  2.1  we  noted  that  the  analytic  current  density  to  a  surface  at  potential  (j)  is  given  by 


j  =  q 


"  2Ee2 

f  E  1 

nr 

J 

f  (E  ±  4>)dE  =  q 


E  ±  (j) 


F(E  ±  (|>)dE 


(61) 


where  L  is  0  for  the  repelled  species  and  l(j)|  for  the  attracted  species,  the  “+”  sign  is  for  ions  and 
the  sign  is  for  electrons,  f(E)  is  the  distribution  function  at  infinity,  and  F(E)  is  the  flux  at 
infinity. 

Below  we  describe  how  these  integrals  are  done. 
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3.1.1  Maxwellian  Currents 


The  Maxwellian  integrals  are  simple. 

|F(E±cj))dE  =  e  ^ 6 


j  =  q 


f  E  A 


VE  ±  4>y 


2710m 


nexp 


E  | 
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—  exp 
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v  9y. 

0  1 

07 

dE 


which  for  the  repelled  species  simplifies  to 


J  =  nei 


and  for  the  attracted  species  simplifies  to 


e9  ^ 
exp 


2?im 


+ - 
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J  =  ne, 


e0 


27im 


V  oy 


v0  y 


The  average  yield  requires  a  numeric  integral. 

Y(E)Eexp 


E  ±  4> 


dE 


Y  = 


Eexp 


E  ±  (j) 


dE 


For  the  attracted  species  use  the  substitution  u  =  exp 


E  ±  (j) 


v  e  j 

-0}  Y(-91nu  +  (|))(-01nu  +  (j))du  }  Y(- 01nu  +  (j>)(-01nu  +  (j>)du 


-0}  (- 01nu  +  (|>)du 


0  +  (j) 


For  the  repelled  species  use  the  substitution  u  =  exp 

^+(j)> 


"_E" 

l  ey 


0exp  ^  }  Y(-01nu)(-01nu)du  J  Y(-01nu)(-01nu)du 

v  y  )i  o 


0exp  —  }  (-91nu)du 

v  0  y> 
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(62) 


(63) 


(64) 


(65) 
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(67) 
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The  numerator  is  calculated  numerically  using  200  points  and  Simpson’s  integration. 

3.1.2  Fontheim  Currents 

The  Fontheim  distribution  includes  both  a  Gaussian  and  a  power  law  term.  This  distribution  is 
only  applicable  for  electrons. 

Gaussian  Term 

The  Gaussian  integrals  are  slightly  more  complex  than  the  Maxwellian. 


j  =  q 


f  E  A 

E-f 


F(E  -  4>)dE  =  q 


f  E  A 

E-<\> 


tC(e  -  <|))exp 


e-^-e/ 


^  v  §auss  J  J 
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(68) 


Use  the  substitution  u  = 


E  4>  E0 


A 


gauss 


j=<q 


C  I 

(AgaUSSu  +  <|)  +  E0)exp(-u2^gaussdu  =JrCqA2gauss  uexp(-u2)du  +  7£q(<|>  +  E0)Agau,syerfc(L)  (69) 


where  L  is  now  max 


'-E„r4>-EA 

gauss  gauss  J 


The  substitution  x  =  exp(-  u2)  can  be  used  to  solve  the  remaining  integral. 


j  =  <qA2gauss  ~exp(-L2)+  <  q((j)  +  E0  )  Agauss  -y  erfc(L) 


(70) 


The  average  yield  requires  a  numeric  integral.  Using  the  same  substitution  and  definition  of  L  we 
have 


oo 

J  Y(AgauSSU  +  <t>  +  E0  XAgaussU  +  <1>  +  Eo  )eXP  (~  U  '  )AgauSsdU 


Y)  =  — 


Agauss  \ exp  (-  L2 )+  (<t>  +  E0  )Agauss  *  erfc(L) 


(71) 


The  integral  in  the  numerator  can  be  divided  into  three  in  order  to  isolate  the  peak  at  u  =  0.  If 
L  >  0,  then  the  first  two  terms  are  zero  and  =  L. 
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(72) 


The  first  and  third  terms  can  be  evaluated  using  the  substitution  x  =  exp(-  u2).  Then 

u  =  +V—  In  x  ,  where  the  sign  is  used  in  the  first  term  and  the  “+”  sign  in  the  third.  Under  the 
assumption  that  £,  is  small,  symmetry  can  be  used  to  write  the  second  term  in  terms  of  the  error 
function.  The  numerator  can  be  written  as 
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exp(-52) 
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exp(-L2) 


Y(-  A gauss +  <t>  +  E„  ]  “  Agauss  +  4=^  dX  +  AgaussY(<l>  +  Eo  +  Eo  ft) 

V  v  —  In  x  y 


A, 


expl 
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V  v  —  In  x  y 


(73) 


Power  Law 


j  =  q 
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(74) 


Using  the  substitutions  x  =  E-(j)  and  then  z  =x'("  1 1  and  z  =x'“  we  get 


j  =  q< 


max(Eu-(J),Eu )  ^ 

dz 


-  -  (a  - 1) 

max  (-<j), El)  ^ 
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+q<<t>  f— 
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q< 
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(75) 


The  same  technique  is  used  to  evaluate  the  average  yield. 
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3.1.3  Kappa  Distribution  Currents 
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which  for  the  repelled  species  simplifies  to 
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for  k  >  1. 

For  the  attracted  species  the  integral  requires  a  substitution. 
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The  average  yield  requires  a  numeric  integral. 
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Use  the  substitution  u  =  1  + 


E±(j) 

kE. 


|y(kEq(u  -l)+  <|))(kE0(u  -l)+  <())(u) ~K_1  du 

(Y)  =  ^ - - - 

J(kEo(U-1)+^XU)  ^du 

I  ' 


(81) 


where  L'  =  1  ±  — for  the  repelled  species  and  L'  —  1  for  the  attracted  species. 
kE., 


Then  use  the  substitutions  y  =  uK  l  and  y=  uK  for  an  even  distribution  of  points. 
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where  L,  = 
the  attracted  species. 

3.1.4  Convected  Maxwellian  Currents 


(82) 


for  the  repelled  species  and  Li  =  L2  =  1  for 


As  the  convected  Maxwellian  distribution  depends  on  the  velocity  of  the  plasma  with  respect  to 
the  spacecraft  and  the  angle  between  the  bulk  plasma  velocity  and  the  particle  velocity,  the  full 
angular  integral  is  needed. 


j,  =  j*  (n,  O,  U)  =  qJJJ  d3  v(  v  •  n)f  s  ( v) 

The  coordinate  system  is  shown  in  Figure  8. 


(83) 
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Figure  8.  Coordinate  System  Used  for  Convected  Maxwellian  Environment,  in  Spacecraft 

Frame 


Written  as  an  integral  over  the  velocity  at  infinity,  v®,  we  have 


Js  =  ne. 


e0 
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2  27t  ti/2 


jd(pjd9jdvoovoo(vl±2^ 
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v*  -2  v^U  cosx+U 


sin0cos0  ,  (84) 


where  the  lower  velocity  limit  L=0  for  an  attractive  potential  and  y[2§  otherwise,  and  %  is  the 
angle  between  the  flow  vector  U  and  the  incident  velocity  at  infinity  v®  that  is  dynamically 
mapped  to  the  incident  angle  0  and  is  given  by 

cos  x  =  cos  0^  cos  0U  +  sin  0^  sin  0U  cos  tp .  (85) 


where  0U  is  the  polar  angle  of  the  flow  velocity  vector  and  0®  is  the  polar  angle  of  v®.  With  this 
substitution,  the  azimuthal  integrand  has  the  lornT  exp(-a  cos(cp)),  and  noting  that  the  zeroth- 
order  modified  Bessel  function  of  the  first  kind  Io  is  given  by 


2  K 

27tl0(a)  =  |  d(pe-acos(p 
0 


(86) 


we  can  write  the  current  as 


oo  n/2 

js  =  Ajdv^v,^  +2c())e‘lV"0  Jd0ev”Ucosa”cos9"Io(vooUsin0oo  sin0u)sin0cos0 

L  0 


with  A  =  qsns 


(87) 


The  velocity  and  polar  angle  integrals  are  performed  numerically.  The  angle  and  magnitude  of 
the  flow  velocity  vector  (defined  in  the  spacecraft  frame)  are  specified  by  the  user,  and  the  angle 


*  The  evenness  of  the  (even-order)  Bessel  function  swallows  a  minus  sign  here. 
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at  infinity  0aj  is  obtained  from  orbit  relations  in  a  1/r  potential,  as  a  function  of  the  total  energy 

V2  s 

E  =  — and  launch  angle  0.  From  reference  1  we  obtain8 


r  2(E  +  c|))sin 2  a 


1+  ll  +  4E^E  +  f^sin2  a  cos  (if  -  9„ ) 


(88) 


Measuring  orbital  angles  from  the  polar  axis  that  is  normal  to  the  sphere  at  the  launch  point  and 
setting  a  =  r  determines  the  angle  90 


cosC^  -90)  =  -f 
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Setting  r  =  oo  gives 
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In  the  limit  as  U— >0  we  recover  the  expected  isotropic/uniform  flux  density,  while  as  the  Mach 
number  increases  the  flux  approaches  a  cosine  dependence  on  relative  particle  angle.  [30] 

For  full  consistency,  the  secondary  electron  emission  current  due  to  ion  or  electron  impact  and 
the  backscattered  current  due  to  electron  impact  would  be  given  by 


00  ,  Til  2 

d9Y,sec,back(v,  e^ev„Ucos3„cos3„  ^  sin  ^  )sin  ^COS  9 

(91) 

As  this  is  a  time  consuming  calculation  and  electron  Mach  numbers  are  always  under  1,  we 
compute  the  yields  due  to  incident  electrons  using  an  isotropic  Maxwellian  distribution. 

•sec, back  /w  \  sec, back  .incident 

Je  \  A  maxwellian /e  Je  (92) 

The  secondary  emission  due  to  ions  is  computed  using  the  full  angular  integral. 

The  effective  density  of  each  component  is  the  expectation  value  of  the  distribution  function.  The 
techniques  used  to  compute  the  average  values  of  the  square  root  of  the  inverse  of  the  energy  are 
the  same  as  used  to  evaluate  the  current  integrals  above. 


Js  = 


Ajdv^v^v2  ±2^^v”  J 


'  Starting  from  Goldstein  eqn  3-46,  substitute  angular  momentum  l  =  rxp  =  rmv  sin(9),  potential  energy  O  =  k/r, 
and  l/2mv2  =  E  +  then  for  a  launch  radius  r  =  a). 
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3.2  Implicit  Charging  Using  the  Boundary  Element  Method 

Among  the  difficulties  of  developing  accurate  and  robust  algorithms  for  spacecraft  charging  has 
been  the  inability  to  calculate  electric  fields  accurately,  let  alone  to  predict  how  electric  fields 
will  change  as  a  result  of  surface  potential  changes.  Nascap-2k  uses  the  Boundary  Element 
Method  (BEM)  [31]  to  calculate  accurate  electric  field  and  as  the  basis  for  implicit  charging 
equations. 


3.2.1  Boundary  Element  Method  Algorithm 

The  Boundary  Element  Method  is  a  means  for  relating  fields  and  potentials  in  a  region  to  sources 
on  the  boundary  of  the  region.  It  is  comparable  to  a  sum  over  the  coulomb  field  of  all  the  charges 
in  a  region  rather  than  an  iterative  field  solution.  In  our  case,  the  region  is  the  space  exterior  to  a 
spacecraft,  and  the  boundary  is  the  spacecraft  surface.  Also,  we  assume  the  “free  space  Green’s 
function”,  i.e.,  the  potentials  in  the  region  obey  Laplace’s  equation. 

The  sources  are  sheets  of  charge  coincident  with  the  spacecraft  model’s  surface  elements.  We 
assume  that  each  surface  element,  j,  has  a  constant  charge  density,  csy  The  familiar  relation  for 
the  potential  of  a  point  charge  then  generalizes  to  an  integral  over  the  object  surface: 


=  - >(4ne0  =Y 

4;rs„r 


d2r.  — 
ru 


(93) 


where  (j);,  which  could  be  the  potential  at  any  point  in  space,  is  considered  the  potential  at  the 
center  of  a  surface  element.  Similarly,  the  familiar  relation  for  the  electric  field  of  a  point  charge 
generalizes  to: 


E  =  -3 - >(te„)E,=y  d!ry(r,-r,) 

4718  „r  r 


(94) 


where  E;  is  the  electric  field  at  some  point  in  space,  the  point  again  taken  at  the  center  of  a 
surface  element. 

We  express  the  relations  of  potential  and  electric  field  to  charge  density  as  matrices: 

*=[^‘1,0,  (E-n),  =Fijtrj  (95) 


These  can  be  combined  to  obtain  a  relation  between  normal  electric  field  and  voltage: 

(E-n),  =FilCkJVJ 


(96) 


This  last  relation  is  the  key  to  developing  relations  between  surface  charge,  surface  potential,  and 
surface  currents  in  order  to  derive  charging  equations. 
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3.2.2  Doing  the  Integrals 

To  get  C'1  we  need  to  do  integrals  of  the  form 


r 


j 


h 


(97) 


There  are  tricks  to  doing  these  integrals,  which  can  be  found  in  the  literature.  Denote  the  field 
point,  r,  as  P,  and  take  the  domain  of  integration  as  triangle  ABC.  Then,  the  vector  from  the  field 
point  to  anywhere  on  the  triangle  can  be  parameterized  as 

iy  =  PA  +  uAB  +  uvBC  (98) 


where  PA,  AB,  and  BC  are  vectors  between  pairs  of  points,  and  u  and  v  are  parameters,  each  in 
the  interval  [0,1].  The  square  of  the  distance,  ry,  then  becomes 


r2  =  r 
y  y 


•rij=  ZTk/u 
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(99) 


0<k,Z<2 


where  the  coefficients  Ty  are  formed  from  pairwise  scalar  products  of  the  three  vectors  above. 
Now,  the  integral  can  be  expressed  as 
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where  the 
The  inner 


Vi  coefficients  are  functions  of  v. 

integral  (over  u)  may  be  found  in  standard  integral  tables  [32] 
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(101) 


We  are  left  to  do  the  outer  integral  (over  v)  numerically.  To  facilitate  this,  we  select  the  vertex  A 
such  that  the  scalar  product  PA-BC  is  the  minimum  of  the  three  choices.  Then,  very  few 
integration  points  are  needed  in  the  outer  integral  (over  v).  (We  use  a  5-point  Simpson’s  rule  in 
the  current  implementation.) 

The  integrals  for  the  electric  field  are  similar,  although  there  are  more  of  them.  The  same 
techniques  apply. 
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3.2.3  Test  for  Accuracy 


We  tested  the  algorithm  for  accuracy  by  calculating  the  electric  fields  on  the  surface  elements  of 
a  uniform  potential  sphere.  The  sphere  model  in  Figure  9  was  originally  defined  using  Patran, 
and  has  90  surface  elements  and  92  nodes.  It  provides  a  fairly  coarse  mesh  over  the  surface  of 
the  sphere. 


Figure  9.  Sphere  Defined  Using  Patran 

The  next  figure  shows  the  result  of  the  BEM  calculation.  Note  that  the  total  variation  of  electric 
field  over  the  surface  is  only  three  percent.  Also,  the  apparent  radius  of  the  sphere  (given  by  V/E) 
is  within  three  percent  of  the  nominal  radius. 
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Figure  10.  Normal  Component  of  Electric  Field  at  Surface  of  Sphere  as  Computed  Using 

the  BEM  Approach 

By  contrast,  the  same  calculation  using  Nascap-2k ’s  finite  element  algorithm  with  higher  order 
elements  gives  an  electric  field  variation  of  seven  percent,  and  a  value  of  electric  field  that  is 
about  fifteen  percent  high.  A  finite  element  algorithm  with  linear  elements  gives  an  eight  percent 
variation,  and  a  value  that  is  about  twenty-five  percent  high. 

3.2.4  Implicit  Charging 

To  develop  charging  equations  we  need  to  express  physical  charges  on  physical  surfaces  in  terms 
of  the  voltages  on  the  same  objects.  We  can  then  proceed  to  calculate  what  voltage  changes  will 
be  produced  by  changes  in  charge  (currents).  Because  the  interior  of  a  spacecraft  is  NOT  free 
space,  the  physical  charge  densities  bear  no  relation  to  the  a’s,  which  we  have  now  eliminated 
from  the  calculation. 

To  get  the  physical  charges  we  use  Maxwell’s  divergence  equation  in  the  form  a  =  V-D.  Figure 
1 1  shows  the  “Gaussian  pillboxes”  used  to  calculate  the  actual  surface  charges  on  insulating 
surfaces  and  on  conductors. 
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Figure  11.  “Gaussian  Pillboxes”  Used  to  Calculate  the  Actual  Surface  Charges  on 
Insulating  Surfaces  and  on  Conductors 


For  an  insulating  surface  (upper  pillbox  in  Figure  11),  the  external  field  is  E-n,  which  we  obtain 
from  the  BEM  solution.  The  internal  field  is  related  to  the  capacitance  between  the  insulating 
surface  and  its  underlying  conductor.  Thus,  the  total  charge,  q,  on  such  a  surface  is  given  by: 


q.  =  Ai(E-n),  +  Cic (4>i  ~<h) 


(102) 


where  A  is  the  surface  area.  For  a  conductor,  since  charges  are  mobile,  it  is  not  useful  to  know 
the  charge  on  each  individual  surface,  but  we  need  to  work  with  the  total  charge  on  the 
conductor.  Conducting  surfaces  include  the  obvious  “bare”  surfaces  (lower  pillbox  in  Figure  11), 
as  well  as  the  surfaces  that  underlie  the  insulating  surfaces  (middle  pillbox  in  Figure  11).  In  both 
cases,  we  have  zero  electric  field  internal  to  the  metal.  To  obtain  the  total  charge  on  all  the 
surfaces  of  the  conductor,  we  need  to  sum  the  external-field  charge  terms  for  the  bare  conducting 
surfaces,  plus  the  capacitive-charge  terms  for  the  insulator-metal  interfaces: 

Qc=y>,(E.n),- 

bare  insulators  nnh 


We  previously  found  the  BEM  expression  for  the  external  fields  in  terms  of  the  cell  potentials: 


(E-n),  =1)^ 


(104) 


Substituting  this  into  the  equations  for  qi  and  Qc,  performing  the  indicated  sums,  and  adding  the 
capacitive  terms,  we  get  the  matrix  equation: 


Q=GO  (105) 

where  the  vectors  are  composed  of  contributions  from  all  the  insulating  surfaces  and  a  single 
term  representing  all  the  bare  conducting  surfaces. 

Q=  {qi,  q2,  ...  ,q„,  Qc}  (106) 

®  =  {<()i,  <|>2,  ((V,  <|>c}  (107) 

where  Q  and  <D  only  contain  entries  for  those  entities  physically  capable  of  accumulating  charge, 
viz.,  conductors  and  insulating  surfaces,  and  the  charges  and  potentials  are  related  by  the 
charging  matrix,  G. 
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We  are  looking  for  an  equation  that  relates  currents  to  voltage  changes.  So,  we  differentiate  the 
charge  equation  in  time: 


I  =  Q  =  GO 


(108) 


Discretize  to  a  finite  time  interval: 


IAt  =  G[0(t  +  At)-0(t)] 


(109) 


Implicitize  by  evaluating  current  at  the  final  time  (also,  for  simplicity,  changing  [t+At,t]  to  [t,0]): 

l(t>  =  G[«(t)-«(0)]  (110 

Linearize  currents  with  respect  to  voltage  (since  we  do  not  know  the  final  voltages  at  which  the 
current  is  to  be  evaluated): 


i(t)*i(o)+r[®(t)-®(o)] 


dm 


And,  solve  the  equation: 


[G  -  I't]o(t)  =  [G  -  I't]o(0)  +  l(0)t  (112) 

This  gives  us  a  straightforward  matrix  equation.  Everything  on  the  right-hand-side  is  known,  and 
we  can  solve  using  standard  linear  algebra  equation  solver  packages. 

Before  proceeding  to  examples,  it  is  worth  commenting  on  the  derivative  of  current, 

i;,=<yaiv 

Consider  three  cases: 


1.  For  current  sources  such  as  incident  plasma  current,  we  usually  approximate  the  current 
as  a  function  of  the  local  voltage  only.  This  gives  a  diagonal  term  in  I'. .  Since  such 

currents  decay  exponentially,  some  care  must  be  taken  not  to  underestimate  the  change  in 
current  if  the  voltage  is  changing  in  such  a  direction  that  a  component  of  current  is 
increasing  (i.e.,  electron  current  for  a  surface  whose  potential  is  increasing  (toward  zero) 
from  a  large  negative  potential). 


For  example,  suppose  we  have  <[)  =  I(4>)  =  1  -  e  <!l  and  cj)(0)=-2.  The  implicit  solution  for 

<K0)(i-rt)+i(<|>(0))t 


<j>(t)  is  4>(  t)  = 


- .  If  we  use  the  actual  derivative!'  =  -exp(cf>(o))we 


(1-I't) 

obtain  the  upper  curve  in  Figure  12,  which  works  satisfactorily  for  sensible  timesteps  like 
0.1  or  1,  but  overshoots  the  solution  (cJ)(qo)  =  0)  badly  for  long  timesteps  like  10  or  1000. 
However,  if  we  recognize  that  zero  is  a  value  we  wish  not  to  overshoot  and  accordingly 
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set  I'  =  -we  obtain  the  lower  curve  in  Figure  12,  which  gives 

0-<|)(0) 

results  for  short  timesteps  and  does  not  overshoot  the  solution,  even  for  a 


equally  good 
timestep  of 


1000. 


0.1  1  10  100  1000 

Timestep 


Figure  12.  Example  Problem  Solutions  for  Short  and  Long  Timesteps  for  Two  Choices  of 

Current  Derivative 


2.  Conductivity  current,  such  as  I;  =  cr(<()c-<j>i),  contributes  off-diagonal  as  well  as  diagonal 


terms  to  the  current  derivative  matrix. 


SL 


SL 

— -  =  -a;  — -  =  +a 


.  Surface  conductivity 


contributes  many  more  off-diagonal  terms. 

3.  The  case  that  is  particularly  problematic  is  I;  =  Ii(E-n).  This  occurs  for  a  photo-emitting 
surface  at  negative  potentials.  The  problem  is  that  the  local  electric  field  is  a  function  of 
the  potentials  on  all  the  surfaces.  However,  the  BEM  provides  us  with  exactly  that 
function.  We  can  now  compute  the  term 


I'  =  dljdty-  =  SI1/5E1  x  cEj/d^ 


(113) 


using  the  relation 


(E-n),  =FlkCk,®J 


(114) 


derived  from  the  BEM. 


In  Nascap-2k  we  impose  some  additional  conditions  on  the  current  derivative  in  order  to  promote 
stability,  some  of  which  are  related  to  the  sharp  cutoffs  that  can  occur  with  secondary  and 
photoemission  currents: 

1.  If  a  surface  is  at  positive  potential  and  has  positive  current  before  any  limiting  of  low 
energy  emitted  currents,  then  the  derivative  is  set  so  as  to  fully  turn  on  the  secondary  and 
photoemission  currents  with  a  potential  drop  of  no  more  than  2.5  V. 
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2.  If  a  surface  is  negative  but  has  positive  current  and  normal  electric  field,  the  current  is 
turned  off  after  rising  to  a  potential  of  no  more  than  +5  V. 

3.  The  current  derivative  is  set  to  allow  a  potential  change  of  no  more  than  2  kV. 

4.  The  current  derivative  is  set  to  allow  a  potential  change  of  at  least  1  V. 

3.2.5  Example:  Sunlit  Sphere 

To  illustrate  the  implicit  charging  model,  we  calculate  the  charging  of  a  sunlit  Teflon  sphere  in  a 
1  cm'3,  20  keV  plasma.  The  NASCAP/GEO  version  of  this  was  published  in  1978.  [33]  A  version 
of  the  original  result  is  shown  below  in  Figure  13.  The  dark  side  of  the  sphere  gradually  charges 
negative  due  to  incident  plasma  electrons,  while  photoemission  grounds  the  sunlit  side.  The  sun 
direction  is  (1,1,1)  (from  the  upper  right).  Eventually,  however,  the  strong  negative  potentials  due 
to  the  dark  surfaces  wrap  around  the  sphere  and  form  a  barrier  to  photoelectron  escape.  The 
potential  of  each  sunlit  surface  is  subsequently  determined  by  the  condition  that  its  electron  field 
has  a  small,  positive  value  (too  small  to  be  seen  in  Figure  13),  so  that  just  the  right  fraction  of 
photoelectrons  can  escape  over  the  barrier. 
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Figure  13.  Potentials  About  Sunlit  Qsphere  in  Space  as  Computed  by  NASCAP/GEO 


Figure  14  and  Figure  15  show  the  BEM  solution  for  a  PATRAN  sphere  and  a  QSphere 
respectively.  The  subsolar  point  is  the  least  negative  point. 
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Figure  14.  Potentials  on  Sunlit  PATRAN  Sphere  in  Space  Viewed  from  Direction  (1,2,3) 
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Figure  15.  Potentials  on  Sunlit  Qsphere  in  Space  Viewed  from  Direction  (1,2,3) 

Figure  16  is  another  view  of  the  BEM  solution,  showing  more  of  the  dark  side.  Despite  the  fairly 
coarse  gridding  on  the  sphere,  the  gradual  potential  gradient  on  the  sunlit  side  and  the  constant, 
large  negative  potential  on  the  dark  side  are  clearly  seen  in  the  BEM  solution. 
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Figure  16.  Potentials  on  Sunlit  PATRAN  Sphere  in  Space 


3.3  Finite  Element  Method  to  Compute  Potentials  in  Space 

In  calculating  the  potential  in  three  dimensions  around  an  arbitrary  object,  Nascap-2k  employs  a 
finite  element  approach  using  right  parallelepiped  elements  and  nonlinear  edge  interpolants. 
Arbitrarily  nested  subdivision  (up  to  seven  levels)  is  used  to  resolve  important  object  features 
while  including  a  large  amount  of  space  around  the  spacecraft  for  determining  accurate  particle 
orbits  and  keeping  memory  requirements  reasonable.  Figure  17  shows  an  example  of  a 
computational  grid  with  progressively  finer  subdivision  to  give  fine  resolution  near  an 
experiment. 
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Figure  17.  Computational  Grid  with  Subdivision 

Nascap-2k  solves  Poisson’s  equation 


-V2(0p/so 


by  solving  the  associated  variational  principle 


0  =  — | 

dV 

i(v<|>)2  +  0 

SHJ 

L2  £oJJ 

c|)V(j)  •  dS 


(115) 


The  first  term  in  the  volume  integrand  corresponds  to  the  Laplacian  operator;  the  second  term  is 
the  space  charge  contribution.  The  surface  contribution,  given  by  the  last  term,  is  used  to  obtain 
the  surface  electric  field. 

Poisson’s  equation  is  solved  with  surface  potentials  fixed.  To  determine  surface  field  values  for 
fixed  potential  elements,  we  use  the  fact  that  the  surface  element  electric  field  is  conjugate  to  its 
potential,  i.e.,  when  the  full  conjugate  gradient  matrix  multiplication  is  performed,  the  surface 
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element  electric  field  (corresponding  to  the  existing  potentials  and  fields)  is  the  residual 
corresponding  to  the  surface  element  potential.  This  is  an  improvement  over  capacitance  matrix 
methods  in  that  it  takes  full  account  of  the  nonlinear  space  charge  in  the  external  space.  In  the 
variational  calculation,  we  use  locally  defined  basis  sets,  that  is,  interpolants  within  each  cubic 
volume  element.  Fine  mesh  volumes  are  given  the  correct  variational  weight,  ensuring  the 
maintenance  of  accuracy  through  mesh  transition  regions.  The  potential  and  electric  field  are 
defined  at  each  grid  node.  The  potential  inside  each  volume  element  is  interpolated  from  the 
values  at  each  of  its  vertices. 


The  volume  integrals  are  solved  volume  element  by  volume  element.  The  potential  (and  electric 
field)  at  each  point  in  the  element  is  the  sum  of  interpolants  times  the  values  at  the  nodes. 


dV  =  -Y 

2. 


(vtfdv,  +•£ 


P*dV. 


(116) 


The  potential  inside  each  element  is  given  by 

4>e(x,y,z)  =  ^Ni(x,y,z)^ 

1  (117) 

where  “i”  are  the  nodes  of  the  element  indexed  by  “e”,  and  the  N)  are  the  interpolants.  Note  that 
each  element  has  four  interpolants  per  node,  whose  coefficients  are  the  potential  and  three 
components  of  electric  field  at  the  node,  as  discussed  in  Section  3.3.1. 

We  then  have 


Vf 


(x,y,z)  =  ^VNi  (x, 


y>z 


Hi 


(118) 


which  allows  us  to  write 


f(vtfdV 


(X ,  XXVNd>6y.dvMx.y.2)MiiidV, 


e  Ve  ia  jp 


(no) 

Z  e  ia  jP 


Where  the  indices  (a,  (3)  run  over  the  potential  and  three  electric  field  components  at  the  nodes 
and  their  corresponding  interpolants.  Note  that  the  Wiajp  depend  only  on  the  shape  of  the  volume 
element. 


The  second  term  can  be  written 


P4» 


J  L£oJ 


=  2  I  -dV  =2  f  X^N„(x,y.zH„dVe  =2>i  Q 


e  J 
V, 


e  J  ia  ^  o 
V„ 


ia  e0 


(120) 


The  minimization  can  then  be  written  as  the  matrix  equation 
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(121) 


_5_ 

64> 


)<pw<i>+-<p  =° 

5^  2  si 


This  function  is  then  expressed  in  terms  of  the  vector,  <p,  of  all  nodal  values  and  a  sparse, 
symmetric,  positive  definite  matrix  W.  The  conjugate  gradient  solution  to  this  problem  requires 
that  we  be  able  to  form  the  matrix-vector  product  W(p,  which  we  can  do  by  adding  up  the 
contribution  to  W<p  from  each  finite  element  without  ever  forming  (or  storing)  the  sparse  but 
extremely  large  matrix  W. 

The  resulting  system  of  linear  equations  is  solved  using  the  Conjugate  Gradient  approach. 

3.3.1  Continuous  Field  Interpolants 

Nascap-2k  uses  a  finite  element  interpolation  scheme  that  has  strictly  continuous  electric  fields. 
This  approach  insures  the  continuity  of  the  electric  field  used  to  compute  particle  trajectories. 
Triquadratic  interpolants  were  also  considered,  but  experience  showed  that  the  difference  in 
weight  between  the  comer  nodes  and  the  edge-center  nodes  leads  to  a  poorly  conditioned  matrix. 

The  basic  interpolation  functions  consist  of  functions  that  have  unit  value  and  zero  slope  at  their 
home  nodes,  and  zero  value  and  slope  at  opposite  nodes: 

F  (x-l)1 

°  l-2x+2x2 


1  °  1  O  ,  O  2 

l-2x  +  2x  (122) 

and  functions  that  have  zero  value  and  unit  slope  at  their  home  nodes  and  zero  value  and  slope  at 
opposite  nodes: 


G„  =  x(l  -  x)F0 
G,  =  x(x  -ljF, 


(123) 


The  coefficient  of  F  corresponds  to  the  potential  at  the  home  node,  and  the  coefficient  of  G 
corresponds  to  the  (negative  of  the)  electric  field  at  the  home  node.  (In  the  element  interior,  both 
F  and  G  are  needed  to  calculate  potential  and  potential  gradient.) 

The  functions  Fo  and  Go  are  shown  in  Figure  18.  It  can  be  verified  that  these  functions  can 
exactly  reproduce  constant,  linear,  and  quadratic  potentials. 
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Figure  18.  Interpolation  Functions  F0  and  G0 

We  generalize  to  three  dimensions  by  assigning  to  each  corner  of  the  unit  cube  four  interpolation 
functions  corresponding  to  the  potential,  x-component  of  potential  gradient,  y-component  of 
potential  gradient,  and  z-component  of  potential  gradient  for  that  node.  The  interpolants 
corresponding  to  node  ijk  are  then  given  by 

N[jk(^y5z)=Fi(x)Fj(y)Fk(z) 

N!jk(x,y,z)  =  G1(x)FJ(y)Fk(z) 

N,"k(x’y’z)=Fi(x)GJ(y)Fk(z) 

Njk(x,y,z)  =  Fi(x)Fj(y)Gk(z)  (124« 

0  12  3 

The  coefficients  of  the  N  interpolants  are  the  node  potentials,  and  the  coefficients  of  the  N 
interpolants  are  the  negative  Cartesian  components  of  the  nodal  electric  fields.  We  elect  to  omit 
the  higher  order  terms  of  the  form  GFG  or  GGG. 


The  potential  inside  each  element  is  given  by 

^(x,y,z)  =  ^N“k(x,y,z>“k 

ijka 


and  the  components  of  the  electric  field  are  given  by 

E(x,y,z)=X^kVN“k(x,y,z) 

ijka 


(125) 


(126) 


3.3.2  Interpolating  Potentials  and  Fields  for  Special  Elements 

Special  elements  are  those  elements  which  contain  surfaces.  Because  surface-containing 
elements  are  not  empty  cubes,  the  potential  interpolation  functions  described  in  the  previous 
section  cannot  be  straightforwardly  applied.  To  develop  finite  element  matrices  for  these 
elements,  we  require  a  prescription  for  calculating  the  potential  and  electric  field  in  the  element 
interior.  Note  that  a  “special  element”  is  bounded  by  three  types  of  surfaces:  (1)  square  surfaces 


45 


bounded  by  grid  edges  and  shared  with  adjacent  (presumably  empty)  elements;  (2)  object 
surfaces;  and  (3)  surfaces  bounded  by  both  object  points  or  edges  and  grid  points  or  edges.  On 
type  (1)  surfaces  we  use  the  potential  and  field  interpolation  described  in  the  previous  section. 

On  object  surfaces,  the  potential  and  electric  field  must  be  expressed  in  terms  of  the  object’s 
surface  element  potentials  and  normal  fields.  Type  (3)  surfaces  must  smoothly  blend  between  the 
two.  We  describe  below  an  algorithm  to  express  the  potential  and  field  at  any  point  in  the 
element  volume  in  terms  of  the  grid  point  potentials  and  fields,  and  the  surface  element 
potentials  and  normal  fields.  The  algorithm  has  the  property  that,  when  applied  to  an  element 
with  six  type  (1)  surfaces,  it  reproduces  the  empty  cube  result  exactly. 

Let  R  be  a  point  in  a  volume  bounded  by  a  set  of  surfaces  { S } ,  each  of  which  may  be  a  triangle 
or  a  planar  convex  quadrilateral.  For  a  given  S,  let  P  be  the  point  on  S  nearest  R,  let  cj)(P)  and 
E(P)  be  the  potential  and  electric  field  at  P,  and  n  be  the  unit  normal  (from  the  surface  point  P 
into  the  volume)  to  the  surface.  Let 


d  =  R-P 
d2  =d*d 

Ks  =  As/d2  ford*n>0 

n=v,  ,k„ 

Z^{s}  s  (127) 


where  As  is  the  area  of  surface  S  and  Ks  is  its  weighting  function. 

Let  /  be  the  distance  from  P  in  direction  n  to  the  next  surface  intersection.  (In  practice,  it  is 
adequate  to  extend  /  to  the  intersection  with  the  surface  of  the  rectangular  parallelepiped 
bounding  the  volume.)  Then  the  contribution  of  S  to  the  potential  at  R  is 

(j)s(R)  =  (Ks  /  N  )[4>(  P  )-(l-d»n//)( d«n )(E(P)  »n)]  (128) 

and  the  total  potential  is  found  by  summing  the  contributions  from  all  surfaces. 

The  electric  field  is  found  by  differentiating  the  above.  The  contribution  of  S  to  the  electric  field 
is 


Es  =  4>s  [VN /  N  -  VKS ) /  Ks )]  +  (Ks  /  N)[-  V<\>(P)  +  (1  -  2d  .  n / l)n(E(P)  •  n)]  (ng) 

Where  V(()( P )  is  taken  tangential  to  the  surface,  and  is  calculated  using  the  continuous  field 
interpolants  of  S  is  a  full  square  in  empty  space,  and  otherwise  using  linear  interpolation  for 
triangles  and  bilinear  interpolation  for  quadrilaterals. 

This  approach  allows  us  to  interpolate  surface  fields  and  potentials  on  bounding  surfaces  in  a 
manner  that  maintains  the  continuous  electric  field  property.  Potentials  and  normal  field 
components  are  defined  at  the  centroids  of  the  original  object  surfaces  and  area-weighted 
averaged  at  surface  comers.  Potentials,  normal  fields,  and  tangential  field  components  at  all 
surface  points  that  are  vertices  of  bounding  surfaces  are  expressed  as  linear  combinations  of 
centroid  potentials  and  fields. 
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3.3.3  Boundary  Conditions 


Most  often  we  use  a  zero  potential  boundary  condition  at  the  outer  boundary  of  the 
computational  grid.  For  some  problems,  specifically,  laplacian  potential  and  linearly  Debye 
screened  potentials,  a  different  approach  is  desirable.  For  these  problems,  we  apply  boundary 
conditions  by  extending  the  problem  beyond  the  computational  box  using  external  elements 
whose  inward  boundary  is  a  quarter-boundary-surface-element  (QBSE)  containing  one  problem 
node  and  which  extends  radially  outward  from  the  center  of  the  computational  box  (to  which  we 
henceforth  refer  as  “the  origin”).  (See  Figure  19)  These  elements  are  added  to  the  sum  overall 
volume  elements.  We  ignore  potential  variation  over  the  QBSE,  assuming  it  to  have  the  potential, 
<(>node,  of  its  problem  node,  so  that  this  element  makes  a  diagonal  contribution  to  W.  Within  that 
approximation  the  new  external  element  subtends  solid  angle  Q.='Ax/R\  where  the  QBSE  is 
normal  to  the  x-axis  at  distance  x  from  the  origin,  and  R  is  the  distance  of  the  center  of  the  QBSE 
from  the  origin.  (R  and  x  are  measured  in  outer  grid  units.)  The  potential  in  the  element  is 


<Kr)  =  <t> 


node 


( 

U  J 


exp  (-(r-R  )//.„), 


(130) 


for  Debye  length  An-  We  can  now  proceed  to  evaluate  the  needed  integrals.  For  the  laplacian 
(A,d=oo)  case  we  get 


i(Vc^)2dV 


/•CO 

r2dr 

f  R^ 
~2 

J  R 

J 

(131) 


so  that  after  properly  taking  account  of  the  mesh  spacing,  L,  the  external  element’s  contribution 
to  the  matrix  element  is 


Wvv  =  )L(x/R2) 

where  x  and  R  are  measured  in  units  of  L. 


(132) 


Figure  19.  An  Extended  Element  with  Analytically  Extrapolated  Potential  Is  Appended  to 
Each  Boundary  Quarter-Square  to  Implement  Non-Zero  (Monopole)  Potential  Boundary 

Conditions 
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For  the  case  of  finite  Debye  length  evaluation  contains  more  terms  and  involves  integrals  that 
cannot  be  expressed  in  closed  form,  so  we  approximate 


dr 

J  R 


e 


-r/  A, 


r 


n 


to  eventually  get 


Wvv  =iL(x/R2)  ^  +  2pln 


2py 


+  ■ 


1 


(L)3x 


1  +  2p  8pA' 


■  + 


(133) 


(134) 


where  p =R  L/X. 

It  is  recommended  that  zero  potential  boundary  conditions  be  used  if  X  is  comparable  to  or 
shorter  than  an  outer  grid  mesh  unit,  and  for  most  particle  tracking  problems  for  which  no  clear 
value  exists  for  the  Debye  length. 

3.4  Space  Charge  Stabilized  Poisson  Iteration 

Poisson’s  equation  can  be  written  dimensionlessly  as 

-V!4>  =  (ni -n  J/A! 


where 


.  (|>  ,  0  A; 

O  =  - ,  and  s  - -  =  — 

0  °  N„eL2 


Debye 


(136) 


A  is  the  dimensionless  Debye  length,  N0  is  the  ambient  plasma  density,  ni  =  Ni/N0,  ne  =  Ne/N0, 
and  the  Laplacian  is  also  normalized  by  L2. 

The  traditional  approach  to  the  solution  of  equation  (132)  has  been  an  explicit  iteration  of  the 
form 


-V2Ov  =  A_2(niOv_1  -neOv_1)  (137) 

where  v  is  the  iteration  index  and  the  charge  density  is  determined  using  the  potentials  of  the 
previous  iteration.  This  method  can  be  shown  to  be  unstable'  when  the  Debye  length,  Ad. 
becomes  small  with  respect  to  other  scale  lengths  of  the  problem.  This  can  be  understood  by 
considering  that  a  smooth  potential  variation  over  a  distance  of,  say,  1000,  would  require  a 

smooth  V2(j)  which  is,  in  turn,  given  everywhere  by  the  charge  density.  Maintaining  a  smooth 
charge  density  distribution  is  difficult  when  any  errors  in  determining  (ne  -  rip  are  multiplied  by  a 
huge  A'  .  There  is  one  effective  remedy  to  this  dilemma  described  in  Parker  [34],  but  the  process 
reported  here  appears  to  be  more  efficient  in  the  short  Debye  length  limit.  This  method  involved 
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the  combination  of  two  concepts.  One  uses  a  partial  implicitization  of  the  repelled  density.  [35] 
The  other  simply  reduces  the  charge  density  to  an  acceptable  level  whenever  the  first  method  is 
inadequate.  An  “acceptable  level”  means  a  value  not  so  great  that  it  will  cause  spatial  oscillations 
in  the  potential. 

3.4.1  Implicitization 

The  right  hand  side  of  Equation  133  is  the  dimensionless  charge  density. 

q(r.®-(r))  =  L-!(n1(®')-ne(®'))  (B8) 

The  charge  density  at  the  present  iteration  may  be  linearized  about  the  previous  potential 
iteration 


q  ( Ov  )  =  q  ( Ov  1 )  +  q'  ( cbv  1 )  (d>v  -  Ov“1 ) 


where  q' 


da 

— - ,  and  the  r  dependence  has  been  dropped  for  clarity.  With  this  expression,  we  may 
00 


write  the  implicit  Poisson  iteration  scheme 


-V2<Dv  - q'  (ov_1 )  <DV  =  q((Dv_1 )  - q'  (®  v“‘ ) ®  V_1 


(139) 


Though  it  is  not  immediately  obvious,  the  implicit  character  of  Equation  135  makes  it  more 
stable  than  scheme  (133)  provided  q'  is  constrained  to  be  non-negative.  This  can  be  understood 
by  realizing  that  in  Equation  133  the  charge  density  is  treated  as  an  independent  variable, 
whereas  in  Equation  135  the  charge  density  is  determined  simultaneously  with  the  potential. 


The  finite  element  approximation  to  Equation  135  produces  the  matrix  equation 


^(we  -S'eVe)®v  =S-S'eOv_1 


(140) 


where  (e)  refers  to  each  element  and  S  is  derived  from  q  by  the  following  analysis. 

3.4.2  Charge  Limiting 

For  A>  1,  S  is  simply  the  total  charge  associated  with  each  node,  q.  (Physically,  this  means  that 
an  entire  grid  cube  of  plasma  with  one  species  eliminated  alters  the  potential  by  no  more  than  the 
temperature.)  However,  for  L«  1,  numerical  noise  and  features  like  a  sheath  edge,  which  may 
span  only  a  few  Xo,  become  incorrectly  amplified  when  the  q  determined  at  a  point  becomes 
multiplied  by  the  entire  nodal  volume.  When  it  is  not  possible  to  reduce  the  volume  element  size, 
stability  can  be  preserved  by  replacing  Q  (and  Q')  with  a  reduced  value  S,  (and  S')  which  is 
calculated  to  be  the  maximum  allowable  charge  for  the  element. 
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Because  of  the  artificial  amplification  argument,  S  is  often  the  more  realistic  total  for  an  element, 
in  the  sense  that  it  produces  potential  variation  appropriate  to  the  spatial  resolution  rather  than 
causing  unphysical  overshoots.  Before  deriving  S,  we  define  the  barometric  potential  Ob  =  ln(ni), 
which  (in  cases  where  electron  density  is  governed  by  the  Boltzmann  factor,  e°)  is  the  potential 
for  which  Q  =  0.  For  the  “Linear”  and  “Non-linear”  space  charge  density  formulations,  the  ion 
charge  density  is  also  dependent  on  the  potential.  For  these  cases,  the  barometric  potential  as 
used  here  is  zero,  Ob  =  0.  Note  that  it  is  important  that  S— >  Q  as  O  — >  Ob  if  quasineutral  regions 
are  to  be  modeled  correctly.  To  determine  S,  consider  a  capacitor  with  potential  difference  Ob  - 
O.  area  L2,  and  a  separation  of  L.  The  charge  qc  on  this  capacitor  is  given  by 


qc  =  cav  =  -p— (0b  —  o)ee 

Le 


(141) 


In  the  units  of  our  previous  q,  qc  becomes 


qLiImt  =a(c£>b  _<i>) 


(142) 


which  is  the  maximum  allowable  charge  per  element,  with  the  parameter  a,  set  to  the  maximum 
value  consistent  with  the  stability  of  the  Poisson  solver.  Thus  at  each  node,  we  choose  for  the  charge 


with 


|S|  =  min  (|qLimit  | ,  |q|) 


(143) 


for  S  =  qLimit 
for  S  =  q 


(144) 


The  way  in  which  S  and  S'  are  computed  for  each  charge  density  formulation  appears  in  Section 
3.4.5  below. 


The  effect  of  this  algorithm  is  as  follows:  If  a  problem  has  been  specified  where  a  boundary 
potential  would  be  screened  in  less  than  an  element  or  two  (the  limit  of  any  code’s  resolution), 
sufficient  sheath  charge  is  redistributed  to  allow  the  potential  to  be  screened  over  the  minimum 
number  of  elements  that  is  consistent  with  stability. 

3.4.3  Analysis  of  the  Space  Charge  Stabilized  Poisson  Method 

The  charge  stabilized  Poisson  method  calculates  for  each  node  the  maximum  allowable  charge 
that  is  consistent  with  the  stability  of  a  linearly  interpolating  Poisson  solver.  This  method  is 
developed  above,  but  a  further  analysis  is  presented  here  to  help  the  user  interpret  its  impact. 

Nascap-2k’s  charge  stabilization  is  accomplished  through  the  process  of  charge  limiting, 
illustrated  in  Figure  20.  This  figure  shows  two  charge  versus  potential  curves  for  the  case  where 
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the  ion  density  is  fixed  (tracked  or  equal  to  the  background  density)  and  the  electron  density  is 
barometric.  Equation  135  is  rewritten  as 

q  =  a  exp  (-<Dm  )  (exp  <f>b  -  exp  <t>)  (M. 

where  Om  =  In (a/.2r)ebye/ L2 )  a  is  user  specified  and  ®b  is  the  barometric  potential,  ®b  =  ln(ni). 
For  curve  qi,  ®b  =  -3.0  (n;  =  0.05);  for  curve  q2,  ®b  =  0.0  (n;  =  1.0)  and  for  both  curves  ®m  =  - 
2.2,  and  a  =  1.  For  each  curve,  the  limiting  charge  as  given  by  Equation  138  is  also  shown.  The 
limiting  charge  is  rewritten  here  as 


q Limit  =  a  (Ob  -O) 


(146) 


which  intersects  the  “natural”  charge  curves  at  Oc  and  Ob.  The  charge  stabilization  method 
reduces  the  charge  to  the  limiting  value  when  O  >  Oc  and  uses  the  natural  charge  for  O  <  Oc. 

The  parameter  Om  provides  a  good  measure  of  the  limiting  process.  From  Figure  20  it  can  be 
seen  that  Om  is  the  point  at  which  the  slope  of  the  natural  charge  curve  equals  that  of  the  limiting 
charge  line.  Figure  21  shows  a  family  of  curves  giving  the  dependence  of  the  cutoff  potential  Oc 
on  the  barometric  potential  for  various  values  of  cpm.  These  curves  were  obtained  by  numerically 
solving  for  the  zeros  of  the  difference  between  q  and  q Limit.  This  difference  equation  always  has 
two  solutions,  one  at  ®c  and  one  at  cpb,  with  the  exception  of  a  degeneracy  at  cpc  =  cpb  =  cpm, 
which  is  indicated  in  the  figure.  This  figure  shows  that  the  charge  limiting  is  minimal  for  Om  >  -1, 
and  quite  severe  for  ®m  <  -6  or  so. 
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Figure  20.  Plots  of  Space  Charge  (Curves  Qi  and  Q2)  as  a  Function  of  Potential  as  Given  by 
Equation  135.  The  Straight  Lines  Represent  the  Maximum  Allowable  Charge  for  Non- 
Oscillatory  Potentials.  The  “Natural”  Space  Charge,  Qi  or  Q2  Is  Acceptable  for  Which 
Slopes  of  the  Curves  and  the  Corresponding  Line  Are  Equal. 


Figure  21.  Plot  of  the  Space  Charge  Cutoff  Potential,  <1>C,  Versus  Barometric  Potential  (<bh  = 
Ln  NO  for  a  Series  of  <bm  Values  (-02,  -0.5,  -1.0,  -2.0,  -3.0,  -4.0  ...  -11.0)  -  the  Point  at 

Which  Om  =  <pb  =  <pc  Is  Also  Indicated 
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3.4.4  Sheath  Boundary  Potential 

Space  charge  limiting  must  be  accounted  for  when  choosing  the  sheath  boundary  potential. 
Space  charge  limiting  limits  the  rate  at  which  space  charge  can  cause  the  potential  to  drop. 

Consider  the  first  volume  element  of  a  sheath  to  satisfy  the  laws  of  Child  and  Langmuir  (planar 
space  charge  limiting).  At  the  sheath  edge  (z=0)  and  one  element  in  (z  =  L)  the  potential  and 
electric  field  are  given  in  Table  1. 

Table  1.  Potential  and  Electric  Field  Variation  Given  by  Planar  Space  Charge  Limiting 


Position 

0 

L 

Potential 

0 

K(L)4/3 

Electric  field 

0 

(4/3)K(L)1/3 

By  Gauss’s  law,  the  charge  per  unit  area  in  this  element  is  given  by 


Q  =  T  K(L)I/3 


s„  A  3 


(147) 


Nascap-2k  computes  the  charge  density  to  be  limited  by  (a/L~)  times  the  mean  potential 
(assuming  linear  interpolation),  and  so  gets 


_Q 

£  A 


=  L 


^  a  ^ 


V ^  J 


(Ml  r) 


(148) 


Equating  these  two  expression  gives  a  =  8/3. 

Because  of  the  requirements  of  disk  space  and  computational  time  of  three-dimensional  codes, 
Nascap-2k  is  frequently  operated  at  high  Om  values,  i.e.,  a  coarse  mesh  with  respect  to  the  Debye 
length.  In  these  cases  charge  is  removed  at  almost  all  points.  Ideally,  the  charge  that  was 
removed  was  excess  charge  generated  by  the  coarse  gridding.  This  is  the  artificial  charge 
amplification  argument.  However,  since  Nascap-2k  must  be  reliably  stable,  the  result  is  that  too 
much  charge  is  often  removed.  This  results  in  an  enlarged  sheath  thickness  for  high  negative  Om 
problems.  In  order  to  compensate  for  this  sheath  enlargement,  the  sheath  boundary  potential  must 
be  adjusted. 

Equating  the  above  equation  to  the  code  resolution,  L,  gives 

fx  =5.1xlCT6L4/3e1/3n2/3  =O.74(L/Xd)4/30  m4c» 
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The  potential  (j)A  may  be  interpreted  as  the  potential  below  which  Nascap-2k  underestimates 
screening.  At  best,  beyond  the  cj)A  contour,  the  potential  drops  about  one  order  of  magnitude  per 
element.  For  9  =  0. 1  eV,  n  =  10  m"  ,  and  L  =  0.2  m,  we  find  <j)A  =  6  V.  If  the  6  V  contour  is 
correctly  placed,  the  0.6  V  contour  lies  at  least  one  element  beyond  (at  the  approximate  sheath 
location),  and  the  default  sheath  contour  (0.07  V)  is  yet  another  element  farther.  This  would 
produce  a  sheath  area  that  is  too  large.  The  suggested  criterion  for  the  sheath  boundary  potential  is 

^sb  =  Max(0  In  2,0.24^)  (150 

where  0.24  =  e"l0/0'7  is  the  planar  screening  per  element  allowed  by  Nascap-2k.  Note,  however, 
that  (j)A  depends  strongly  on  the  grid  in  which  the  sheath  is  found,  so  that  if  an  increase  in  object 
potential  moves  the  sheath  from  grid  3  to  grid  2,  a  corresponding  larger  value  of  c()sb  should  be  used. 

3.4.5  Charge  Density  and  Derivative  in  Nascap-2k 

Nascap-2k  has  a  selection  of  eight  different  expressions  to  describe  the  charge  density  for 
different  kinds  of  problems.  Each  expression  uses  implicitization  and  charge  limiting  to  a 
different  extent.  As  there  is  no  space  charge  for  the  Laplace  space  charge  option,  the  following 
discussion  does  not  apply  to  that  case.  The  expressions  used  for  the  space  charge  and  its 
derivative  including  limiting  are  listed  below. 

The  parameter  D  is  the  local  mesh  spacing,  L,  divided  by  the  user  provided  Debye  limiting 
value,  which  defaults  to  2.  The  parameter  g  is  the  maximum  of  plasma  density  reduction  factor 
computed  by  neutral  wake  model  (0<g<  1 )  and  10  6.  The  basic  strategy  in  the  formulations  below 

<3P  £  (!) 

is  that  Debye  length  analog.  A,  defined  by  — —  =  — V ,  must  be  no  greater  than  D. 

d(j>  X 


To  see  the  relation  between  L,  D,  and  A,  consider  the  solution  of  the  Helmholtz  equation  in  a 
single,  linear  zone.  The  variational  principle  for  the  Helmholtz  equation  is 


0  = 


5 

5^ 


(V*)2  + 


1*1 

2  A2 


(151) 


The  squared-gradient  term  is  minimized  by  a  constant  potential,  whereas  the  squared-potential 
term  is  minimized  by  going  symmetrically  through  zero.  So,  there  must  be  some  minimum  value 
of  A  for  which  the  solution  does  not  change  sign.  Minimizing  the  variational  function  on  the 
interval  [0, 1]  with  respect  to  4>(1)  when  4>(x)  =  (1-x)  4>(0)+x  4>(1)  gives 


<Ki)  =  <|>(0) 


2  — 


3A2 


2  +  - 


3A 


(152) 
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so  that  the  condition  for  non-oscillatory  potentials  becomes  X2  > 


\J_ 

~6 


where  we  have 


redimensionalized  to  the  local  mesh  spacing,  L.  The  Debye  limiting  value  following  from  this 
derivation  would  be  V6  =  2.45  . 


Laplace.  The  Laplacian  space  charge  option  specifies  that  the  charge  density  is  zero  and  no 
space  charge  limiting  is  needed. 


(153) 

i.e.,  charge  exists  only  on  object  surfaces  and  external  boundaries,  as  determined  by  the 
boundary  conditions.  Space  charge  iterations  may  still  be  required,  however,  due  to  the  treatment 
of  surface  electric  fields. 


d(p/£o)_0 

(154) 

Linear  (Debye  Shielding).  The  Linear  space  charge  option  solves  the  Helmholtz  or  Debye  - 
Huckel  equation: 


P_  =  _A 

So  ^nl 

^nl2=  max(?4/g,D2) 

d(pAQ_  1 


(155) 


Non-Linear.  Non-linear  space  charge  is  appropriate  for  most  low-earth-orbit  type  plasmas. 
Poisson’s  equation  is  solved  with  space  charge  given  by: 


p/e0=-(^/X2) 


max(l,C((j),E)) 


i+^/eJ3 

C((|),E)  =  min((Rsh/r)2 ,3.545 


(Rsh/r)z=2.29|E^nl/0nl|1'262|0nl/^ 


|3/2 

nl| 

0.509 


Xnl2=  max()^p  /g,D2) 

0nl= 


(156) 


The  term  C(<|)E)  (analytic  focusing)  comes  from  fitting  a  finite  temperature  spherical  (Langmuir- 
Blodgett)  sheath.  If  analytic  convergence  is  turned  off  (on  the  Potentials  Advanced  screen)  then 
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C(4>,E)  is  set  to  zero.  Additional  adjustments  are  made  to  the  convergence  portion  of  this  formula 
for  significant  values  of  velocity. 


d(p/e0) 

d(j) 


=  Max 


p/s0 


max(l,C((|),E)) 


1  +  M^4tt  4>/0n 


|3/2 


v2  J  1  +  ^/4^  |(j)/0n]|3/' 


(157) 


Note  that  the  derivative  of  the  convergence  contribution  to  the  charge  density  with  respect  to  the 
potential  is  assumed  to  be  negligible. 


Frozen  Ion.  The  Frozen  Ion  formulation  is  intended  for  short  timescale  (typically  sub¬ 
microsecond)  problems  for  which  it  is  a  good  approximation  to  assume  that  ions  remain 
stationary  and  at  ambient  density  (“ion  matrix”  approximation),  but  electrons  achieve  barometric 
equilibrium.  The  space  charge  function  depends  on  the  mesh-dependent  potential,  cjq  <  0,  which 
satisfies 


l-expf^/e)  = — (a,/d)2  (<h/e) 


(158) 


for  D>a.  Otherwise,  cj)i=-lxl0'6.  The  meaning  of  this  quantity  is  that  the  charge  density  increases 
linearly  at  a  stable  rate  for  negative  potentials  between  zero  and  cjq,  and  then  asymptotes  to  the 
ambient  ion  density  for  potentials  more  negative  than  (jq. 

The  space  charge  is  then  given  by 


p/s0  = 


(^i/E>2)(1-exp(^i))  ^>0 

-^/D2  0>(|)>^ 

-^i/D2+(0/^2ebye)(exp(^i/0)-exp(V0)) 


(159) 


d(p/s0) 

d(|) 


Max 

l/D2 

Max 


p/£0 


,ex 


p(4>/4>i)/d2 


"^fL’exp(^,/e)/^°cbye 


<\>>o 

0  >  4>  >  4>, 


(160) 


The  derivative  of  charge  with  respect  to  potential  appears  in  the  variational  function  as  a 
coefficient  of  cj)  .  If  this  is  too  large,  then  minimizing  (j)  means  you  oscillate  from  one  grid  point 
to  the  next  in  order  to  get  lots  of  intermediate  zeroes,  which  is  not  a  satisfactory  solution. 

Full  Trajectory  Ions.  Ion  densities  are  calculated  from  steady-state  ion  trajectories.  Electrons 
are  barometric. 
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(161) 


—  =  —  (l  -  exp  (min  (10,  (<|>  -  c|>b  )/0ft ))) 
£0  £0 


pf=max(p““,en„n) 

<^b  =01nl 


9ft  =  max 


V  en  y 
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V  So  J 
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d(p/go) 

dcj) 


e0  <>  — <l>b 


Pi 

8.9, 


max 


(<t>-4»b)/0ft  > 1 
abs((<|)  —  <|)b  )/0ft  )  <  10 
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p  1 

pf  exp(min(lO,(^-^b)/0ft))A 

So^-^b’ 

So0|t 

otherwise 


(162) 


Hybrid  PIC.  This  algorithm  is  used  for  timescales  (typically  sub-millisecond)  on  which  it  is 
practical  to  treat  ion  motion,  but  electrons  may  be  considered  in  barometric  equilibrium.  The  ion 
density  is  computed  from  actual  ion  macroparticles.  The  electron  charge  density  is, 
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Full  PIC.  For  this  option,  it  is  assumed  that  the  total  charged  density  (electron  and  ion)  was 
previously  computed  and  stored. 


„  _  „  tracked 

P-P 

d(P/£o)_Q 

dc|) 


(165) 


3.5  Macroparticles 

Macroparticles  represent  either  current  or  charge.  Macroparticles  that  represent  current  are  used 
for  steady-state  problems  and  ones  that  represent  charge  are  used  for  dynamic  problems. 
Macroparticles  are  also  used  for  visualization. 

3.5.1  Macroparticle  Creation 

When  macroparticles  are  created  at  a  sheath  surface,  the  sheath  surface  is  divided  into  triangles 
and  a  macroparticle  representing  the  current  from  that  section  of  the  sheath  is  created  at  the 
center  of  the  triangle. 

When  macroparticles  are  created  at  the  problem  boundary,  the  outer  surface  of  each  grid 
boundary  volume  element  is  divided  into  squares  (of  size  requested  by  the  user)  and  a 
macroparticle  is  created  at  the  center  of  each  square  to  represent  either  the  current  through  that 
area  (“Boundary”  and  “BField”)  or  the  charge  passing  through  that  area  during  a  user  specified 
time  (“Boundary  Injection”). 

When  macroparticles  are  created  throughout  the  volume  (for  uniform  charge  density 
initialization  or  charge  exchange),  they  are  created  at  points  appropriate  to  the  weighting 
functions. 

When  macroparticles  are  created  at  user  specified  locations  with  user  specified  velocities,  each 
macroparticle  represents  current. 

For  visualization  only,  macroparticles  can  also  be  created  at  the  intersection  of  the  sheath  surface 
and  a  user  specified  plane. 

3.5.2  Macroparticle  Tracking 

Macroparticles  are  tracked  using  a  third  order  energy  conserving  algorithm. 

x(t  +  dt)  =  x(t)+  v(t)dt  +  — dt2  +  — dt3  +  — dt4 
v  ’  w  v  ^  2  6  24 

v(t  +  dt)=  v(t)+x,dt  +  —  dt2  +  — dt3 

2  6  (166) 


where 
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(167) 


x2  =  — (E  +  vxB) 
m 

x3  =  —  (v*VE  +  x,  xB) 
m 

x4  =  —  (x3  xB) 
m 

The  electric  field  is  that  at  the  location  of  the  macroparticle  determined  by  interpolation  from  the 
nodes  of  the  volume  element.  The  magnetic  field  is  a  constant,  uniform  field  supplied  by  the 
user. 

Macroparticles  are  tracked  for  a  timestep,  which  can  be  as  long  as  a  minute  for  static  problems  or 
as  short  as  microseconds  for  dynamic  problems.  The  timestep  is  divided  into  substeps  such  that 
the  macroparticle  can  travel  no  more  than  a  fraction  of  the  local  mesh  size  in  a  single  substep.  By 
default  the  fraction  is  0.1;  the  user  may  specify  an  alternative  value.  The  particle’s  total  energy  is 
calculated  at  the  beginning  of  the  timestep  and  preserved  during  the  timestep.  Changes  in 
potential  between  timesteps  will  be  reflected  by  change  in  the  particle’s  total  energy. 

If  the  macroparticle  represents  charge,  at  the  end  of  each  timestep,  and  optionally  at  the  end  of 
each  substep,  the  charge  and  current  are  assigned  to  the  nodes  of  the  element  the  macroparticle  is 
in  using  the  same  weights  as  the  potential.  The  charge  density  (charge  divided  by  element 
volume)  is  assigned  to  the  element.  If  the  assignment  is  performed  at  the  end  of  each  substep,  the 
charge  assigned  is  divided  by  the  fraction  of  the  timestep  that  the  substep  represents. 

If  the  macroparticle  represent  currents,  at  the  end  of  each  substep  its  charge  density  (current 
times  substep  time)  is  optionally  assigned  to  the  element  within  whose  bounds  the  macroparticle 
is. 

When  a  macroparticle  impacts  a  surface  element,  its  current  is  assigned  to  the  surface  element. 

3.5.3  Macroparticle  Splitting 

Macroparticles  are  split  both  to  avoid  having  heavy  macroparticles  in  well-resolved  regions  and 
to  simulate  a  thermal  distribution.  Macroparticle  splitting  has  been  implemented  for 
macroparticles  read  from  an  external  file,  for  a  uniform  distribution  of  macroparticles  for 
initialization,  and  for  macroparticles  injected  from  the  boundary.  Macroparticles  can  also  be  split 
when  entering  a  more  finely  resolved  grid. 

General  Principles 

1.  Macroparticles  are  split  in  velocity  space  only.  Because  we  frequently  find  ourselves  in  high- 
field  regions,  spatial  splitting  would  raise  problems  with  energy  conservation. 

2.  To  be  split  in  velocity  space,  a  macroparticle  must  carry  a  temperature.  We  assume  the 
temperature  is  always  isotropic.  The  fission  products  carry  half  the  temperature  of  the 
original  macroparticle,  while  the  remaining  thermal  energy  appears  as  kinetic  energy  of  the 
split  macroparticles. 
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3.  For  splitting  purposes,  we  define  the  Z-axis  to  be  along  the  direction  of  the  macroparticle 
velocity,  the  X-axis  randomly  chosen  in  the  plane  normal  to  Z,  and  the  Y-axis  mutually 
perpendicular. 

4.  We  split  into  two  or  three  macroparticles  with  added  velocity  along  each  axis,  except  that  we 
may  elect  not  to  split  along  the  Z-direction  if  the  kinetic  energy  exceeds  the  thermal  energy. 
Not  splitting  along  Z  helps  ameliorate  macroparticle  proliferation,  but  makes  an  error  by  not 
preserving  the  original  macroparticle  temperature  along  Z.  We  thus  end  up  with  eight,  nine, 
or  twenty- seven  new  macroparticles. 

5.  Macroparticle  velocity  is  assumed  to  be  acquired  by  acceleration  rather  than  actual  drift  (i.e., 
spacecraft  velocity).  If  there  is  actual  drift  (e.g.,  ram  velocity),  then  the  drift  velocity  should 
be  removed  before  splitting  the  macroparticle,  and  added  back  after. 


6.  If  the  macroparticle  with  speed  uo  is  split  by  two  along  the  X  or  Y  axis,  the  new  velocity  is 
+  0.707A/e6/rn  .  Along  the  Z  axis,  the  velocity  increment  is  calculated  as  if  the  temperature 


were  0-2mu; 


1  +  - 


mu; 


—  1 


7.  If  the  macroparticle  with  speed  uo  is  split  by  three  along  the  X  or  Y  axis,  there  is  a  zero- 
velocity  central  macroparticle  and  two  “probe”  macroparticles  with  velocity  ±O.866A/e0/m  . 
Along  the  Z  axis,  the  velocity  increment  is  calculated  as  if  the  temperature  were 


0  -  2mu ; 


1  +  - 


0 


mu; 


-1 


Examples  of  the  effects  of  macroparticle  splitting  appear  in  Plasma  Interactions  with  Spacecraft, 
Final  Report. 

3.6  Particle-in-Cell  Tools 


3.6.1  Transverse  Surface  Currents 

We  developed  a  pseudopotential  approach  (similar  to  the  velocity  potential  used  to  describe 
potential  flow  [36]  in  fluid  dynamics)  to  the  solution  of  the  current  continuity  equation  for  the 
computation  of  transverse  surface  currents  on  a  spacecraft  acting  as  an  antenna.  As  a  boundary 
condition,  one  surface  element  of  each  antenna  element  is  specified  as  connected  to  the  biasing 
power  supply,  and  the  solution  provides  the  required  current.  The  vector  transverse  surface 
current  in  each  surface  element  is  that  which  provides  the  best  fit  to  the  edge  currents. 

A  solution  having  no  circulating  currents  is  guaranteed  if  we  assign  a  pseudopotential  value  to 
each  element  of  surface  and  take  the  current  across  their  common  edge  to  be  proportional  to  the 
difference  in  their  pseudopotentials.  (In  that  case,  current  that  flows  “downhill”  would  need  to 
flow  back  “uphill”  in  order  to  circulate.)  The  current  is  also  taken  as  proportional  to  the  edge 
length  and  inversely  proportional  to  the  distance  between  the  centers  of  the  elements. 
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The  result  of  the  above  treatment  is  a  matrix  that  multiplies  the  vector  of  pseudopotentials  to 
describe  the  buildup  of  charge  in  the  surface  elements,  which  is  then  set  equal  to  the  source 
vector  generated  from  Nascap-2k  results  at  each  timestep.  This  system  of  equations  is  then 
solved  using  the  ICCG  (Incomplete  Cholesky  Conjugate  Gradient)  [37]  algorithm. 

Figure  22  shows  graphically  the  quantities  that  make  up  the  surface  current  equation  for  each 
surface  element.  Moving  the  surface  currents  to  the  left  side  of  the  equation,  and  all  other 
(known)  quantities  to  the  right,  and  expressing  the  edge  currents  as  proportional  to 
pseudopotential  differences,  results  in  a  sparse  symmetric  matrix  that  is  amenable  to  solution 
using  the  ICCG  algorithm  once  the  surface  element  connected  to  the  power  supply  (denoted  the 
“injection  element”)  is  set  to  a  fixed  potential.  The  current  required  to  bias  the  conductor  may  be 
obtained  by  evaluating  the  equation  associated  with  the  injection  element  and  can  be  verified  to 
be  the  sum  of  the  total  change  in  charge  less  plasma  currents  for  all  the  remaining  surfaces  of  the 
conductor.  Surface  currents  are  solved  separately  for  each  electrically  isolated  component. 


Total  change  in  charge  Plasma  Currents  Surface  Currents 


Figure  22.  The  Change  in  Charge  on  a  Surface  Element  Is  Made  Up  of  Plasma  Currents 

and  Surface  Currents 


3.6.2  Propagating  Fields 

The  transverse  surface  currents,  the  volume  ion  currents  computed  during  particle  tracking,  and 
the  volume  electron  currents  saved  in  the  database  by  an  external  code,  are  a  source  of 
propagating  electromagnetic  fields.  Nascap-2k  can  compute  the  magnetic  field,  the  vector 
potential,  and  the  rate  of  change  of  the  vector  potential  from  these  currents. 


A(r)  =  Z 


(168) 


B(r)  =  E 


jitransx(r-r1)AI 


(r-r,)2 


y  Jk"  x (r  - rk )Vk  ^  jfc  x (r -rk)Vk 
k  (r-rk)2  k  (r-rk)2 


where  the  sums  are  over  the  surface  elements  (i)  and  volume  elements  (k). 


(169) 
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LIST  OF  SYMBOLS 


b; 

Range  parameters 

d 

Thickness  or  distance 

e 

Magnitude  of  the  electron  charge  (1.60  x  10"19  C) 

f 

Distribution  function 

g 

Reduction  in  plasma  density  in  wake  in  the  absence  of  surface  potentials 

i  and  j 

indices 

Jth 

Plasma  thermal  current 

jx 

component  of  the  current 

l 

length 

m 

Particle  mass 

n 

Component  density 

na 

Ambient  plasma  density 

n 

Surface  normal 

q 

Charge  on  a  surface  or  particle  charge 

qi 

Range  parameters 

r 

Distance 

t 

time 

v  or  y 

Velocity 

u  or  v 

parameters  in  BEM 

x  or  x 

Position  or  path  length 

y 

substitute  integrand 

A 

Area  or  Point  of  triangle  in  BEM  discussion 

A 

BEM  matrix 

B  and 

B 

Magnetic  field  and  magnetic  field  magnitude  or  point  of  triangle  in  BEM  discussion 
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c 

Capacitance  or  point  of  triangle  in  BEM  discussion  or  Convergence  factor 

Dcl 

Child  Langmuir  distance 

D 

electric  field  (in  Maxwell  equation) 

E 

Particle  energy  in  eV 

E  and 

E 

Electric  field  and  electric  field  magnitude 

E0 

Parameter  of  Gaussian  component  of  Fontheim  distribution  function 

F 

^max 

Energy  at  which  secondary  yield  due  to  normally  incident  electrons  peaks 

F 

Flux  or  interpolant 

F 

BEM  matrix 

G 

BEM  matrix 

H 

Heaviside  step  function 

G 

interpolant 

I 

Current  (amperes) 

I 

Current  vector 

L 

Mesh  spacing  or  lower  limit  of  energy  or  velocity  integral 

Nj 

Interpolant 

P 

Point  of  triangle  in  BEM  discussion  or  point  in  volume 

Q 

Finite  element  matrix 

R 

Resistance  or  Sphere  radius  or  Escape  depth  (range)  or  point  in  volume 

S 

Surface  area  or  Stopping  power  or  Surface  index  or  reduced  charge  density 

AT 

Time  interval 

T 

pairwise  scalar  products  in  BEM 

U 

Velocity  of  plasma  with  respect  to  spacecraft,  -Vsc 

V 

Volume  or  BEM  coefficients 

Vsc 

Spacecraft  velocity 

w 

Finite  element  matrix 
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Y 

Yield 

Z 

Atomic  number 

a 

Parameter  used  in  ion  secondary  yield  or  Parameters  used  in  space  charge  limiting  or  parameter  used  to 
specify  power  law  component  of  Fontheim  distribution  function 

p 

Parameter  used  in  ion  secondary  yield 

X 

Angle  between  flow  vector  and  the  particle  velocity 

^max 

Peak  yield  for  normal  incidence 

Permittivity  of  vacuum  (8.854  x  10'12  F  m"1) 

8 

Relative  dielectric  constant 

4> 

Surface  potential  or  Space  potential 

4>b 

Barometric  potential 

11 

Backscatter  coefficient  or  Fraction  of  the  distribution 

K 

Parameter  used  in  kappa  distribution  function 

xD 

Debye  length 

V 

iteration  index 

e 

Plasma  temperature  in  eV 

p 

Resistivity  or  Charge  density 

a 

Surface  charge  density 

dCEX 

Charge  exchange  cross  section 

\ 

Small  number  used  to  evaluate  current  integral 

V 

Incident  angle  or  Angle  between  neutralizer  axis  and  look  direction 

C 

Coefficient  in  Fontheim  formula 

^gauss 

Parameter  of  Gaussian  component  of  Fontheim  distribution  function 

A(j) 

Difference  in  potential  across  dielectric  layer 

<D 

Potential  vector 

n 

Solid  angle 
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